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INTRODUCTION 

In  the  last  twenty  years,  the  operating  environment  of  farm  tractors  has 
changed  drastically  in  size  and  ergonomics.  The  environment  has  become  quite 
comfortable,  particularly  on  the  larger  tractors.  Often,  operators  tend  to 
use  their  more  comfortable  tractors  to  perform  even  the  light-duty  tasks  which 
do  not  utilize  the  full  power  of  the  big  mp.  i.-s.  When  a  large  tractor  is 
operating  under  a  light  load,  fuel  can  potentially  be  saved.  Most  operators 
tend  to  operate  the  engine  at  a  higher  engine  HPM  than  necessary;  if  they 
would  gear  up  the  transmission  and  lower  the  RPM,  they  could  move  to  a  more 
fuel-efficient  location  in  the  engine  load  map  while  still  producing  the  same 
power  output. 

To  this  end,  extension  personnel  at  Kansas  State  University  have 
encouraged  an  adjustment  in  operating  practices,  called  gear  up  and  throttle 
back.  A  concern  expressed  by  some  operators  about  this  change  in  driving 
habits  was  that  they  might,  inadvertently,  throttle  the  engine  down  to  a  range 
that  could  damage  it  under  sustained  low  speed  operation.  To  avoid  such  an 
occurrence,  a  gear  selection  aid  was  developed  (Blumanhourst .  1984),  (Blu- 
manhourst,  et  al.  1984).  The  device  monitored  the  power  being  used  and 
displayed  a  new  engine  speed,  gear  ratio,  and  fuel  savings  if  more  than  .5 
gal/hr  of  fuel  could  be  saved  by  switching  from  the  present  to  the  new  set- 
ting. The  device  was  generally  successful  at  saving  fuel,  although  the  amount 
saved  depended  on  the  operating  habits  of  the  person  and  the  nature  of  the 
load.  The  .5  gal/hr  figure  was  used  as  the  minimum  to  prevent  the  algorithm 
from  constantly  recommending  minor  changes  in  gears  and  engine  speeds  for 
insignificant  fuel  savings.  The  use  of  the  threshold  implies  that  potential 
fuel  savings  still  exist  even  if  the  device  recommendations  are  followed 


faithfully.  Blumanhourst  (1984,  Fig.  6)  shows  that  more  fuel  can  be  saved  by 
increasing  shift  frequency.  However,  in  the  range  of  20-30  shifts/hr,  the 
fuel  savings  is  small  compared  to  the  inconvenience  of  shifting  often. 

Histograms  obtained  by  Blumanhourst  showed  variations  both  within  a  field 
operation  and  between  field  operations.  For  some  operations  where  the  load 
varies  only  slightly,  many  of  the  load  changes  may  be  below  the  level  that 
would  translate  to  significant  fuel  savings  if  the  gears  were  shifted.  For 
the  more  abrupt  load  changes  common  to  primary  tillage  tools,  fuel  savings 
could  also  be  difficult  to  realize  due  to  the  impossibility  of  shifting  fast 
enough  to  maximize  the  savings. 

The  present  project  seeks  to  capture  the  potential  fuel  savings  by  using 
a  continuously-variable  transmission  to  more  closely  approximate  the  point  of 
maximum  fuel  efficiency.  The  laboratory  setup  consists  of  a  Caterpillar  3304 
engine  coupled  to  a  Vadetec  continuously-variable  transmission  (CVT)  in  series 
with  a  Funk  power-shift  transmission  and  planetary  gearbox  with  a  Midwest 
eddy-current  dynamometer  loading  the  system.  The  entire  system  is  under  com- 
puter control  so  that  the  engine  will  remain  in  the  most  fuel-efficient  area 
of  the  engine  map. 

To  test  the  control  algorithm  for  stability  and  actual  fuel  savings,  a 
"typical"  pattern  of  loading  due  to  tillage  patterns  will  be  needed.  Data 
files  from  the  gear  selection  aid  project  were  used  to  provide  the  test  pat- 
terns. The  many  hours  of  collected  data  were  to  be  condensed  and  represented 
by  a  short  10-15  minute  sample.  Methods  such  as  selecting  a  short  piece  of 
data,  using  the  increment  from  the  last  reading,  or  deviation  about  a  local 
mean  were  considered.   Finally,  it  was  decided  to  represent  the  pattern  using 


statistical  methods  such  as  spectral  density  or  autoregressive-moving  average 
processes.  Thus,  a  short  cycle  could  be  characteristic  of  the  hours  of  data 
collected  in  the  field. 

OBJECTIVE 

The  objective  of  this  study  was  to  represent  various  tillage  implements 
by  characteristic  power  series  obtained  from  data  of  previous  fieJd  trials. 

LITERATURE  REVIEW 

Standards  and  Soil  Models 

For  guidelines  on  how  to  derive  a  typical  test  cycle  simulating  field 
tillage  on  a  dynamometer,  one  potential  source  is  the  tractor  test  standards 
published  by  the  engineering  societies.  In  the  United  States,  the  American 
Society  of  Agricultural  Engineers  (ASAE  S209.5)  and  the  Society  of  Automotive 
Engineers  (SAE  J708  JUN80)  use  the  same  standard  to  test  the  claims  of  the 
manufacturer.  These  tests  subject  the  tractor  to  drawbar,  PTO,  and  fuel  con- 
sumption tests.  The  testing  is  done  on  a  hard  surface  at  a  constant  speed  for 
consistency  and  convenience.  The  International  Organization  for  Standardiza- 
tion also  lists  standards  (ISO  789/1)  which  are  similar  to  the  ASAE  and  SAE 
standards  for  testing  tractors. 

Relationships  betwleen  the  variables  which  affect  a  tractor's  performance 
on  a  hard  surface  may  relate  test  track  data  to  actual  working  conditions. 
Information,  including  graphs,  on  variables  such  as  the  wheel  loading  factor, 
tractor  size  and  configuration,  and  tractive  performance  on  a  concrete  surface 
is  given  by  Leviticus  and  Reyes  (1985).   In  addition,   ASAE  standard  D230.3 


gives  guidelines  to  predict  field  performance  from  test  data  and  estimate  the 
expected  draft/unit  cross  section  for  various  tools. 

The  above  test  procedures  all  deal  with  a  constant  level  of  power.  These 
methods  work  satisfactorily  for  checking  durability  and  estimating  average 
fuel  consumption.  However,  the  load  histograms  obtained  by  Blumanhourst 
(1984)  show  that  the  load  does  not  remain  constant,  but  varies  as  the  tractor 
progresses  through  the  field. 

The  change  in  draft  force  while  pulling  tools  through  the  field  perhaps 
could  be  calculated  by  using  soil  models  to  give  a  test  cycle  which  imitates 
the  loading  on  a  tool.  However,  these  models  are  still  in  an  elementary  stage 
of  development.  Grisso  and  Perumpral  (1985)  evaluate  different  soil  models, 
all  of  which  use  the  basic  wedge  in  front  of  the  tool  to  calculate  the  draft 
force  and  make  other  simplifying  assumptions  such  as  uniform  soil  which  would 
tend  to  make  the  load  constant. 

Industry  Simulations  from  Data 

Industry  and  government  use  road  simulations  to  test  a  variety  of  vehi- 
cles without  the  expense  of  field  testing.  The  Environmental  Protection 
Agency  prescribes  a  sequence  of  automobile  speeds  to  check  exhaust  gas  emis- 
sions (SAE  J1094a)  using  a  dynamometer.  On  actual  roads,  the  road  profiles 
can  be  examined  (see  Chaka  (1978)  for  a  description  of  an  electronic  road  pro- 
filer) to  verify  that  no  changes  in  input  have  occurred  to  vehicles  on  dura- 
bility tests.  Chaka  also  says  that  the  road  profiles  are  useful  in  duplicat- 
ing public  road  profiles  on  special  test  roads. 


6 

Histograms  have  been  used  to  simulate  the  conditions  of  a  tractor  in  the 
field  (Whelpley,  1973).  Field  data  of  axle  torque  was  collected  and  analyzed 
to  produce  frequency  of  occurrence  histograms  of  various  field  operations.  A 
histogram  using  20  load  increments  was  constructed  for  each  unique  field 
operation.  The  individual  histograms  were  then  combined  on  a  time-weighted 
basis  to  produce  master  histograms  for  a  specific  operation  (i.  e.  third-gear 
operation  using  heavy  draft  implements).  These  histograms  were  then  used  to 
make  master  histograms  for  the  engine  and  the  drive  train.  The  field  data  was 
combined  based  on  percent  of  time  spent  at  heavy,  medium,  and  light  draft 
loads  to  make  a  final  histogram  consisting  of  ten  blocks.  For  the  engine  his- 
togram, data  from  all  gears  was  used  in  proportion  to  the  amount  of  actual 
field  usage.  For  the  power  train,  the  histogram  consisted  only  of  data  for  a 
specific  gear  ratio.  The  load  levels  for  the  dynamometer  were  then  selected 
through  a  randomized  block  sequence.  The  final  histogram  could  also  be 
represented  in  continuous  form  by  the  amplitude  probability  distribution 
(Bekker,  1969),  which  uses  a  single  curve  to  display  the  load  levels.  As  an 
alternative  to  the  histograms,  Whelpley  also  considered  using  a  mean  torque 
signal  combined  with  a  random  signal  to  simulate  the  deviation  of  the  load 
history. 

Cryer  and  Nawrocki  (1976)  used  spectral  density  to  determine  the  inputs 
for  a  road  simulator.  Four  basic  steps  were  used  in  this  approach.  First, 
the  vehicle  was  driven  on  a  road  of  interest  while  instruments  recorded  the 
vertical  acceleration  on  each  wheel  spindle.  Second,  the  vehicle  was  placed 
on  appropriate  actuators  and  the  dynamic  characteristics  which  relate  the 
spindle  response  to  the  actuator  inputs  were  measured.  Third,  using  the 
desired  response  obtained  from  the  test  road,  the  effective  road  inputs  were 


derived  by  using  the  dynamic  characteristics  of  the  test  stand.  Fourth,  the 
response  to  the  derived  inputs  was  determined  and  compared  to  the  desired 
response.  The  actuator  inputs  were  then  varied  iteratively  until  the  power 
spectral  density  plots  for  the  simulated  road  profile  matched  that  obtained 
for  the  actual  surface.  Cryer  and  Nawrocki  (p.  2)  report  that  "the  PSD 
description  of  road  profiles  is  gaining  acceptance,  since  it  is  an  accurate 
description  of  the  road"  and  that  it  has  been  proposed  to  the  ISO  as  a  stan- 
dard description  of  road  roughness. 

Time  Series  Models 

Evenly-Spaced  Methods 

Several  methods  of  analyzing  time  series  are  widely  encountered  in  the 
literature.  Two  closely-linked  methods  are  autocovariance.  commonly  displayed 
as  autocovariance  versus  lag,  and  spectral  analysis  which  maps  the  relative 
covariance,  known  as  power  spectral  density,  versus  frequency.  The  autoco- 
variance shows  patterns,  if  any,  as  the  time  variable  changes  while  power 
spectral  density  shows  whether  the  spectrum  has  peaks  of  variance  in  narrow 
frequency  ranges,  as  would  be  the  case  with  a  cosine  curve,  or  whether  the 
relative  covariance  changes  gradually  over  a  wide  frequency  range.  The  clas- 
sic text  on  spectral  density  analysis,  and  with  it  autocovariance,  was  written 
by  Blackman  and  Tukey  (1958).  The  theory  is  now  covered  in  many  books  that 
deal  with  time  series  analysis.  The  basic  theory  is  outlined  in  Appendix  A 
along  with  a  procedure  used  by  Walls  et  al .  (1954)  and  Wendenborn  (1966)  to 
generate  PSD  plots  for  discrete  data. 

Another  type  of  time  series  analysis  consists  of  the  AutoRegressive  and 
Moving  Average,   ARMA(p,q),   models   which  Box   and  Jenkins  (1976)  deal  with 
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extensively.  These  models  consist  of  three  parts,  the  autoregressive ,  the 
moving  average,  and  the  random  plant  noise.  The  autoregressive  part  makes  the 
present  value  dependent  on  p  past  values,  each  multiplied  by  a  distinct  coef- 
ficient which  reflects  the  degree  of  dependency.  The  moving  average  part 
makes  the  present  value  dependent  on  q  past  random  noises.  The  random  noise 
term,  present  in  all  models  even  if  q  is  zero,  is  added  to  reflect  the  assump- 
tion that  the  next  value  is  not  determined  absolutely,  but  will  fall  within  a 
certain  range,  assumed  to  be  gaussian  in  distribution. 

Irregularly-Spaced  Methods 

The  theory  for  time  series  analysis  methods  like  power  spectral  density 
or  the  ARMA  models  has  been  well-developed  for  continuous  or  evenly-spaced 
data  (see  Otnes  and  Enochson  (1972)  or  Box  and  Jenkins  (1976)).  Only  recently 
have  methods  been  suggested  for  dealing  with  irregularly-spaced  data.  Various 
methods  for  analyzing  irregular  time  series  are  given  by  Parzen  (1984). 

Irregular  spacing  can  be  dealt  with  by  weighting  the  individual  elements. 
Masry  (1984)  uses  a  mean  value  of  sampling  rate  which  is  a  constant  in  his 
spectral  density  calculations.  Marquardt  and  Acuff  (1984)  weight  the  i.jth 
element's  contribution  to  the  spectral  density  by  the  proportion  of  time 
represented  in  a  square  containing  that  element.  An  attempt  was  made  to 
analyze  the  available  data  using  this  method  (see  "Methodology").  The  Direct 
Quadratic  Spectrum  Estimation  (DQSE)  method  developed  by  Marquardt  and  Acuff 
is  summarized  in  Appendix  B. 

Harvey  and  Pierse  (1984)  describe  an  ARMA  method  that  is  useful  when  the 
data  was  taken  at  long  intervals  initially  and  then  later  taken  at  closer 
intervals,  as  would  be  the  case  with  stock  prices  which  were  recorded  yearly 


9 

at  first,  and  later  changed  to  monthly  index.  Their  technique  reverses  the 
series  so  to  obtain  the  initial  covariance  matrix.  The  state  space  approach 
is  then  used  to  evaluate  the  parameters  with  an  optimization  program  which 
uses  numerical  derivatives  to  fit  the  parameters.  They  also  show  an  approach 
for  evaluating  an  AutoRegressive  Integrated  Moving  Average  (ARIMA)  method  for 
use  with  differencing,  the  change  between  the  data  points  instead  of  the 
actual  points,  and  seasonal  differencing. 

Jones  (1980,  1985a)  details  an  ARMA  method  which  is  useful  for  cases  of 
evenly-spaced  data  with  missing  values.  The  method  uses  the  state  space 
approach  on  a  given  data  set.  The  initial  covariance  matrix  is  calculated 
recursively  from  the  expected  data  values,  which  can  be  expressed  in  terms  of 
the  coefficients  of  past  values  and  random  noise.  The  scheme  of  model  fitting 
begins  with  initial  estimates  of  the  parameters  and  then  uses  a  nonlinear 
optimization  program  to  vary  the  parameter  values  until  the  minimum  value  of 
-2  ln( likelihood)  is  reached.  This  method  is  described  further  under  the  sec- 
tion entitled  "Theory  of  the  ARMA  Process".  Jones  also  considers  the  ARMA 
method  for  the  case  of  data  with  irregular  spacing  (Jones.  1985b). 

In  the  ARMA(p,q)  model,  which  bases  the  present  value  upon  p  past  values 
and  q  past  random  plant  noises,  the  accuracy  of  the  estimate  of  the  p  and  q 
coefficients  used  with  the  past  values  can  be  judged  by  any  of  several  guide- 
lines (see  Jones  (1985a)  for  a  summary).  The  least  squares  approach  selects 
the  value  that  minimizes  the  sum  of  squares  of  the  residual.  The  Yule-Walker 
equations  choose  the  parameters  which  best  predict  the  present  covariances 
from  the  past  covariances.  The  maximum  entropy  approach  is  similar  to  least 
squares,  except  that  a  backwards  autoregression  is  considered.  Forward  and 
backward  autoregressions   converge   to   the   same   value.    The   term   in  the 
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numerator  for  the  parameters  is  the  same  for  both  directions,  but  the  denomi- 
nator differs  by  a  term.  Burg's  maximum  entropy  estimate  uses  an  average  of 
the  two  denominators  from  forward  and  backward  autoregression  as  the  denomina- 
tor. 

Yet  another  way  of  estimating  the  best  fit,  under  the  assumption  that  the 
noises  are  normally  distributed,  is  the  likelihood  function  which  depends  upon 
a  combination  of  the  innovation;  that  is,  the  difference  between  the  actual 
and  predicted  values,  and  the  variance  of  the  innovation.  The  likelihood 
function  is  often  simplified  to  ln( likelihood) .  The  relative  value  of  likeli- 
hood for  a  data  set  is  considered  rather  than  the  absolute  since  the  value 
differs  according  to  the  type  and  length  of  the  data  sample.  Akaike  (1974) 
makes  another  modification  to  the  likelihood  function  by  adding  a  penalty  for 
increasing  the  order  of  the  autoregression.  Akaike 's  Information  Criterion 
(AIC)  is  meant  to  be  an  estimate  of  2N  times  the  average  information  for  dis- 
tinguishing between  the  assumed  and  true  model,  where  N  is  the  number  of  data 
points.  The  first  part,  consisting  of  -2  ln( likelihood) ,  is  the  penalty  for 
the  badness  of  fit  between  the  modeled  value  and  the  actual  value  in  the  data. 
The  second  part,  the  addition  of  twice  the  number  of  parameters,  is  the 
penalty  for  increased  unreliability  caused  by  using  an  increasing  number  of 
variables  to  describe  the  system.  In  the  concept  of  the  AIC  number,  the  best 
model  is  the  one  with  the  lowest  AIC  value.  Models  with  an  AIC  value  close  to 
the  lowest  can  also  be  considered,  with  a  lower  order  model  being  preferable 
to  a  more  complicated  model. 
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ORIGINAL  DATA  COLLECTION 

The  raw  data  collected  in  the  original  gear  selection  aid  project  con- 
sisted of  time,  engine  and  transmission  hertz  signal  counts,  and  position  of 
the  injector  pump  rack.  The  transmission  ratio  and  ground  speed  were  thus 
easily  obtained.  The  power  and  fuel  flow  were  obtained  by  regression  equa- 
tions from  rack  position  and  rpm.  The  specific  fuel  consumption  was  then 
easily  calculated.  The  above  values  were  recorded:  the  variables  considered 
in  the  present  analysis  are  time,  velocity,  and  power. 


The  instrumentation  sampled  the  data  points  at  a  rate  of  1  sample/sec. 
At  the  time  of  the  experiment,  the  interest  was  in  the  overall  patterns  so  a 
long  time-length  sample  was  desired  to  accurately  estimate  the  load  histogram. 
Tape  was  chosen  as  a  storage  media,  so  in  order  to  work  within  the  storage 
space  of  the  microcomputer  and  minimize  the  tape-writing  time,  the  sampled 
values  were  recorded  in  the  cycle  described  below  and  shown  in  Fig.  1. 
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Figure  1 .   Data  Cycle 


1.  Readings  were  taken  each  second. 

2.  For  the  first  minute  of  the  5-minute  cycle,  the  data  was  stored  as  15  4- 
second  averages.  The  data  treated  in  this  way  will  be  referred  to  as 
fine  data. 
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3.   The  subsequent  4  minutes  of  readings  were  stored  as  4  1-minute  averages. 
The  data  treated  in  this  manner  will  be  referred  to  as  coarse  data. 

The  short-term  (4-sec)  averages  showed  the  faster  variations  while  the  60-sec 
averages  showed  the  long-terms  variations  in  loading  patterns.  The  result  of 
the  data  was  used  to  plot  histograms  of  power-usage  variation  (Blumanhourst , 
1984).  The  data  files  contained  several  segments  of  data  since  the  original 
project  disregarded  values  when  the  gears  were  shifted  or  the  tractor  was 
traveling  from  one  location  to  another.  The  longest  "continuous"  data  segment 
was  used  to  represent  each  file.  A  FORTRAN  program  checked  the  fine  and 
coarse  values  for  the  beginning  and  end  of  the  segment.  The  end  occurred  if 
values  were  missing  when  the  total  accumulated  time  was  less  than  or  equal  to 
5  minutes  or  when  the  count  was  above  30  minutes  and  the  time  skipped  was 
greater  than  or  equal  to  10  minutes.  The  total  sample  length  varied  by  loca- 
tion and  ranged  from  30  minutes  to  over  4  hours. 

MATERIALS 

The  materials  used  in  the  current  analysis  consisted  of  the  data  files, 
collected  as  described  above,  computers  to  analyze  the  data,  and  programs  to 
control  the  data  processing.  A  DEC  PDP  11/34  was  used  to  run  the  first  two 
methods  of  analysis  using  a  spectral  density  approach.  A  Harris  800  was  used 
to  run  Jones'  exact  likelihood  method  due  to  the  speed  of  the  Harris  and  the 
amount  of  iteration  necessary  in  the  program.  The  three  methods  considered 
for  analysis  are  described  below. 
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METHODOLOGY 

Analysis  of  the  series  was  attempted  using  three  different  methods.  The 
first  two  attempts  were  based  on  the  idea  that  any  pattern  can  be  represented 
as  a  collection  of  sine  and  cosine  waves.  The  third  method  of  analysis 
assumed  that  the  present  value  is  dependent  on  the  past  values  as  well  as  on 
random  noises. 

Regularly-Spaced  Spectral  Density 

The  first  analysis  attempt  used  Blackman  and  Tukey's  power  spectral  den- 
sity which  deals  with  regularly-spaced  data  with  no  missing  values  (Walls,  et 
al,1954)  and  (Wendenborn,  1966).  Since  the  available  data  was  evenly  spaced 
with  missing  values  in  the  coarse  segments,  it  was  forced  into  evenly-spaced 
i:i.  . .:  :  s  by  recombining  the  data.  The  basic  interval  used  was  the  4-second 
interval.  The  mean  of  each  minute  of  data  was  determined.  For  the  first 
minute  of  the  data-taking  cycle,  this  was  the  mean  of  15  points.  For  the 
remaining  four  minutes,  this  was  simply  the  recorded  value,  since  this  was 
already  the  average  over  60  seconds.  The  data  was  then  recombined  starting 
with  the  first  minute  of  fine  data.  The  deviations  from  the  one-minute  mean 
of  the  fine  data  were  determined.  The  data  was  then  reorganized  by  combining 
the  deviations  of  the  fine  data  from  minutes  1,6,11,16,21,...  about  the  means 

of  minutes  1,2,3,4,5 This  recombined  data  was   then  analyzed  by  the 

autocorrelation-fourier  transform  method  described  in  Appendix  A.  Since  there 
was  little,  if  any,  statistical  support  for  the  above  recombination  of  data, 
the  search  was  continued  for  methods  that  were  suitable  for  the  data  spacing 
pattern. 
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Direct  Quadratic  Spectrum  Estimation 

Marquardt  and  Acuff's  Direct  Quadratric  Spectrum  Estimation  (UQSE)  method 
was  used  in  the  second  approach  (Marquardt  and  Acuff,  1984).  This  method  cal- 
culates the  spectrum  by  giving  a  weighting  factor,  which  is  dependent  on  the 
time  since  the  last  sample,  to  each  correlation  value.  The  samples  are 
assumed  to  be  point  samples  taken  at  irregular  intervals.  Since  the  algorithm 
had  the  capability  of  dealing  with  irregular  spacing,  the  samples  were 
analyzed  on  a  distance  basis,  as  determined  from  the  average  speed  at  each 
data  sampling,  instead  of  on  a  time  basis.  The  data  in  the  present  case  is 
not  a  point  sample,  but  rather  an  average.  The  fine  samples  are  close  enough 
to  be  regarded  as  point  samples,  and  indeed,  the  apparent  amplitudes  in  aver- 
aged data  and  the  point  samples  can  be  related  by  regression.  The  minute 
averages  are  too  different  from  the  4-sec  averages  to  be  regarded  as  point 
values  and  so  they  are  left  out  of  the  spectral  density  analysis  by  the  DQSK 
method.   The  Marquardt-Acuf f  algorithm  is  given  in  Appendix  B. 

The  most  obvious  result  of  the  weighting  by  time  intervals  is  that  the 
beginning  and  end  values  of  the  minute  of  fine  data,  within  the  sample,  are 
given  a  much  larger  weight  than  the  other  values  in  the  minute.  Consequently, 
known  curves  in  the  form  A  cos(2*;tx/T)  where  A  is  the  amplitude,  T  is  the 
period,  and  t  is  the  time  were  also  analyzed  by  this  method.  The  known  curve 
was  treated  in  the  same  manner  as  the  original  data.  Values  were  calculated 
at  1  second  intervals.  For  the  first  minute,  15  values  were  recorded,  each 
consisting  of  averages  over  4  seconds.  The  next  four  minutes  were  not 
recorded,  as  they  are  not  representative  of  the  same  area  as  the  4-second 
value.  The  cycle  was  repeated  until  a  suitable  length  was  obtained.  On  the 
spectral  density  plot  of  a  known  curve,  an  abrupt  peak  occurs  at  the  expected 
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frequency.  On  either  side  are  prominent  side  lobes,  their  presence  due  to 
leakaf^e  and  their  height  to  the  effect  of  irregular  spacing.  These  lobes  make 
interpretations  of  the  resulting  spectrum  inexact.  In  order  to  obtain  a  sharp 
window  which  focuses  only  on  the  frequency  under  consideration,  as  recommended 
by  Marquardt  and  Acuff,  the  lag  was  long  enough  to  stretch  over  several  5- 
minute  cycles. 

An  attempt  was  made  to  relate  the  main  peak  in  the  spectral  density  plot 
of  a  known  curve  to  the  original  amplitude  by  regression  equations.  The  fac- 
tors affecting  the  height  were  mainly  the  amplitude,  A.  of  the  original  curve, 
and  to  a  lesser  degree,  the  length  of  the  data  sample  and  the  period,  T,  of 
the  wave  (due  to  the  averaging).  The  regression  equations  showed  a  poor 
degree  of  correlation  with  a  changing  rate  of  variance  so  this  approach  was 
reexamined  by  limiting  the  weight  assigned  to  each  data  value.  With  the  ori- 
ginal Marquardt-Acuf f  approach,  the  weight  was  assigned  to  the  minimum  of  the 
length  of  the  Nyquist  wavelength  and  the  halfway  distance  between  the  data- 
taking  intervals.  This  resulted  in  the  starting  and  ending  points  of  the  data 
being  given  much  greater  weight  than  any  other  values.  Limiting  the  weight 
that  could  be  assigned  to  each  point  by  setting  the  maximum  possible  distance 
to  10  m,  which  was  approximately  the  distance  represented  in  each  4-second 
interval  in  the  data  samples,  eliminated  most  of  the  sharp  peaks  and  gave  a 
spectrum  that  was  quite  flat.  Since  the  influence  of  the  beginning  and  end 
values  of  the  fine  data  samples  could  be  unpredictable,  yet  another  method  was 
sought. 
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Autoregressive-Moving  Average  Approach 

The  third  and  last  means  of  analyzing  the  data  is  the  autoregressive  mov- 
ing average  (ARMA)  approach.  Jones  (1980,1985a)  shows  an  adaptation  of  the 
ARMA  process  to  deal  with  missing  data.  A  way  to  recursively  estimate  the 
parameters  is  (Jones, 1980)  by  changing  to  the  state  space  form  as  derived  by 
Akaike  (1974).  The  fit  of  the  model  is  evaluated  by  minimizing  the  -2 
ln( likelihood)  function.  The  model  calculates  likelihood  where  the  data 
exists  and  ignores  portions  of  nonexistent  values.  Thus,  the  value  is  exact 
for  the  given  data.  See  the  section  "Theory  of  the  ARMA  Process"  for  more 
detail . 

The  model  fitting  was  done  by  using  a  set  of  subroutines  written  by  Jones 
(1978).  The  main  program  was  based  on  one  supplied  by  Jones  and  modified  to 
read  in  the  data,  figure  the  average  and  standard  deviation,  and  subtract  the_ 
average  from  the  data  values  to  yield  numbers  distributed  about  zero.  As  in 
the  DQSE  approach,  only  the  fine  values  were  used;  the  coarse  were  neglected. 
The  program  checked  40  models  in  the  form  ARMA(p,q),  with  varying  integers  as 
p  and  q  values.  The  AR  value  p  ranges  from  0  to  6  and  the  MA  value  q  from  0 
to  5.  The  sum  of  p  and  q  is  ^  10.  due  to  program  limitations,  so  the  highest 
order  model  is  ARMA(6,4).  When  the  order  of  p  or  q  is  increased,  the  value  of 
the  transformed  new  variable  is  set  to  zero  and  the  parameter  values  from  the 
previous  model  are  used  to  initialize  the  present  model.  If  a  problem  is 
encountered  in  finding  the  minimum  -2  In(likelihood)  in  the  previous  model, 
such  as  exceeding  the  iteration  limit,  the  array  of  transformed  parameters  is 
reset  to  zero  for  the  next  model.  The  subroutines  calculate  the  -2 
ln( likelihood)  function  and  then  feed  the  parameters  into  a  non-linear  optimi- 
zation system  which  varies   the  transformed  parameter  values  until  the  -2 
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ln{ likelihood)  is  minimized.  After  the  optimization  routine  has  completed  its 
work,  the  main  program  prints  the  output  and  begins  the  next  model.  After  all 
40  models  have  been  evaluated,  the  valid  solutions  are  sorted  into  ascending 
order  of  AIC  (Akaike's  Information  Criteria  is  the  sum  of  -2  In(llkelihood) 
plus  twice  the  number  of  parameters  used  in  the  model.)  and  the  orders  of  p, 
q,  their  coefficients,  estimated  error  variance,  and  AIC  are  printed. 

The  optimization  program  used  in  this  experiment  was  developed  at  the 
National  Center  for  Atmospheric  Research  in  Boulder,  Colorado.  It  finds  the 
local  minimum  of  a  twice  continuously-dif ferentiable  real-valued  function  f 
from  a  given  starting  point  x  .  The  input  for  this  subroutine  consists  of  the 
dimension,  n,  of  the  problem,  a  subroutine  to  evaluate  the  function.  in  this 
case  the  -2  ln( likelihood)  function,  and  an  estimate  x  of  the  minimum  of  the 
function.  When  the  program  is  finished,  the  values  returned  are  an  approxima- 
tion to  the  local  minimum,  the  value  of  the  gradient,  and  a  flag  to  designate 
the  reason  for  stopping.  Termination  may  not  necessarily  be  due  to  conver- 
gence to  the  solution;  other  reasons  may  be  errors  in  input,  inefficient  use 
of  the  program  by  setting  the  dimension  of  the  program  to  1,  or  exceeding  the 
iteration  limit. 

Some  minor  changes  were  made  in  the  routines  obtained  from  .Jones.  The 
Markrep  routines  called  for  IMSL  routine  LEQTIF  to  be  used  as  a  linear  equa- 
tion solver  using  Gaussian  elimination  (Crout  algorithm)  with  equilibration 
and  partial  pivoting.  The  IMSL  library  was  not  available  so  a  routine  based 
on  the  Cholesky  method  was  used  instead  ( James ,et . al . ,  1977).  The  subroutines 
for  calculating  -2  ln( likelihood)  were  modified  so  that  they  checked  the  sign 
of  the  (1,1)  element  in  the  covariance  array,  which  should  always  be  positive, 
before  attempting  to  take  the  square  root  of  that  element  (dimensioned  real). 
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According  to  Jones,  a  negative  sign  could  be  caused  by  roundoff  error  and 
could  be  alleviated  by  inserting  a  small  observational  error  into  the  reading. 
Consideration  of  observational  error  was  available  in  the  routines,  but  not 
used  in  order  to  contain  the  number  of  models.  The  program  stopped  prema- 
turely on  6  of  24  data  files  before  this  check  was  implemented.  Roundoff 
error  was  also  observed  in  one  test  case  where  the  values  were  small  and  many 
were  missing,  in  a  pattern  of  15  present,  60  missing.  A  solution  was  ade- 
quately obtained  after  doubling  the  values  in  the  file  and  processing.  The 
non-linear  optimization  routine  also  was  modified  to  check  the  sign  of  numbers 
before  taking  the  square  root.  No  case  was  observed  where  this  error 
occurred,  but  the  check  was  inserted  as  a  precaution. 

The  output  of  the  optimization  routine  was  also  modified.  Originally, 
the  program  had  three  levels  of  messages  which  would  print  out  varying  levels 
of  information  as  the  program  progressed,  none  of  which  would  print  out  no 
information  and  also  allow  the  model  to  run  an  ARMA(0,1)  or  ARMA(1,0)  model. 
Since  the  program  could  generate  sizable  amounts  of  monitoring  data  while  run- 
ning 40  models,  the  desired  change  was  made. 

THEORY  OF  THE  ARMA  PROCESS 

The  following  is  the  theoretical  basis  for  the  subroutine  set  Markrep, 
which  calculates  -2  In(likelihood) ,  as  given  by  Jones ( 1980 , 1985a) ,  and  is 
included  here  for  the  benefit  of  the  reader.  An  autoregressive  moving  average 
process  with  an  order  of  p,q  (ARMA(p,q))  and  a  zero  mean  can  be  expressed  as 


^i  =  ,^Vt-k  ^  \  ^  ,^  ^k^-k      •  fiJ 

k=l  k=l 

The  first  summation  term  represents  the  autoregressive  part  where  the  present 
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X  is  dependent  on  p  previous  x's.  The  second  summation  term  determines  the 
contribution  that  the  past  shocks  make  to  the  present  vaJue  of  x;  only  the 
past  q  values  are  considered,  hence  the  model  is  referred  to  as  ARMA(p,q). 
The  e  term  is  a  part  of  the  model  no  matter  what  degree  of  q  is  used.  The 
"real  x"  is  considered  to  be  specified  within  a  certain  range  which  is  indi- 
cated by  the  magnitude  of  6  .  The  random  noise  is  the  only  input  from  the 
environment  into  the  system.  The  value  has  a  defined  mean,  in  this  case  0, 
with  the  true  value  distributed  about  the  mean.  Thus  even  an  ARMA(O.O)  model 
would  not  always  have  a  value  of  0:  the  mean  of  the  data  would  be  0  and  the 
numbers  would  be  distributed  within  a  defined  variance. 

An  assumption  is  made  that  the  process  is  stationary;  in  other  words,  if 
the  process  is  sampled  from  time  t  to  t+k  and  then  sampled  at  time  t+m  to  t+n, 
the  mean  will  be  the  same  and  the  variance  of  the  shock  noise  will  be  the 
same.  The  values  of  k,  m,  and  n  are  arbitrary  and  the  value  of  the  mean  is 
not  necessarily  zero  for  the  definition  of  stationarity  to  be  met.  The 
assumption  of  stationarity  means  that  the  roots  of 

P      k 

1  -   E  a,  z   =  0  [2] 

k=l   " 

must  lie  outside  of  the  unit  circle.   The  coefficients  of  the  x     (equation 

t  —  K 

(1))  not  have  to  lie  between  -1  and  1  to  meet  this  condition.  For  example, 
the  roots  of  the  equation  could  be  2  and  5/4.  both  of  which  lie  outside  of  the 
unit  circle;  equation  (2)  would  then  be 

1  -  1.3z  +  .4z^  =  0 

The  equation  must  also  be  invertible;  that  is,  the  equation  below  must  have 
roots  outside  the  unit  circle. 


1  +  I  /3  z  =  0 
k=l  " 
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[3] 


Without  this  constraint,  the  influence  of  the  random  shocks  would  continue  to 
grow  as  time  progressed  and  the  uniqueness  of  the  parameters  would  be  lost. 

The  state  of   the  process   in   Akaike's   Markovian   representation 
(Akaike, 1974)  is  defined  as  a  column  vector,  Z,  with  a  length  of  m,  where 


m  =  Max(p, q+1 ) 


Z(t) 


X(ti t) 
X(t+l! t) 


[4J 


[5J 


.x(  t+m-1  It). 
where  x(t+jlt)  is  the  prediction  j  steps  ahead  when  given  the  x  values  up  to 

the  time  t.   Hence,  x(t|t)  is  simply  the  value  x    Rewriting  equation  (1)  for 

a  one-step  prediction  using  this  terminology  gives 


x(t.lit)  =   E  Vt.l-k^  ,\    Vt.l-k 
k=l  k=l 

Expanding  to  a  j-step  prediction  uses  the  previous  predictions  and  yields 


[6] 


J-1 


x(t+jlt)  =  Z  aj^x(t  +  j-kit)  +  E  aj^x^^.  j^ 


+  Z   (3,  £ 


,  ^k  t*-j-k 


[7J 


k=l  k=j      ■'  k=j 

If  the  indices  are  beyond  their  proper  range,   the  summations  concerned  are 
eliminated.   Updating  equation  (7)  to  time  t+1 ,  it  becomes 

J-2  p  q 

x(t+j!t^l)  =  I  a  x(t+j-klt  +  l)  +   I   «wX      +   Z   K^t^i   ic  f^J 

k=l  ^  k=j-l  ^  ^^'^       k-^j-1  ^  ^^'^ 

Subtracting  equation  (7)  from  (8)  gives  a  recursion  which   involves  only  the 
random  input  at  time  t+1: 
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J-1 


x(t+jit+l)  -  x(t+jit)  =  L     a,  [x(t+j-kit-l)  -  x(t+j-k't)]  ^   /3 .  ,e^  , 

k=l  ^  J""^  ^^^ 


[91 


If  the  term  in  brackets  is  redefined, 


x(t-j|t+l)  -  x(t+j!t)  =  gjet^.1 


10] 


where  g.  is  some,  as  yet  undefined,  constant.   Equation  (9)  can  be  rewritten 
as 


x(t+jlt+l)  =  x(t+jit)  + 


J-1 


J-1 


k=l 


k^j-k 


't  +  1 


[11] 


The  quantity  in  brackets  may  be  combined  into  one  general  variable  called  g, 
where 


^1  =  ' 


[12] 


and  /3  =  0  for  j  >  q.   Equation  (11)  can  then  be  written  more  compactly  as 


x(t+jit  +  l)  =  x(t^jlt)  +  g-e^.^j 


[14- 


with  the  final  element  in  the  state  vector  written  as 


x(t+mit+l)  =  Z  a,  x(t+m-k!t)  +  g  e  , 
,  ,   k  m  t-^1 

k=l 


[15] 


The  autoregressive-moving  average  (ARMA)  model  can  be  described  in  a 
recursive  fashion  by  using  the  state  space  framework.  As  an  example,  the 
technique  is  shown  for  a  first-order  autoregression  and  later  generalized  to  a 
higher-order  model.   The  equation  for  the  ARMA(1,0)  model  consists  of: 
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x(t+l)  =  ax(t)  +  e(t)  [16J 

where  x(t)  is  the  state  at  time  t,  assuming  zero  mean,  a  is  a  vaiue  represen- 
tative of  the  state  transition  matrix  (For  this  first-order  process  tho  state 

transition  matrix  is  a  scalar.),  and  e(t)    is  a  random  shock  with  a  variance  of 

2 
a    .      If  the  process  is  observed  with  error,  the  observation  equation  becomes 

y(t)  =  x(t)  +  v(t)  [11] 

where  x(t)  is  the  "true"  process,  v(t)  is  a  white  noise  process  with  a  vari- 
ance R,  and  y(t)  is  the  process  for  which  values  are  actually  obtained.  This 
modification  of  (16)  is  equivalent  to  an  ARMA(1,1)  process  since  y(t+l)  now 
depends  on  the  previous  value  at  (t)  and  an  error  in  addition  to  a  random 
excitation.   The  variance  of  this  estimated  state  is 

P(t+lit)  =  E{[x(t+1)  -  x(t+llt)]^}  [18] 

At  the  beginning  of  a  zero-mean  process,  the  best  estimate  of  the  state 
is  zero: 

x(0!0)  =  0  [19] 

The  variance  of  the  process  at  the  beginning  is  best  estimated  as 

2 

P(0!0)  =  — ^— T  [20] 

1  -  a 

This  is  the  variance  that  would  be  expected  from  a  no-skills  forecast.  For 
the  special  case  of  first-order  autoregression,  the  recursion  proceeds  as  fol- 
lows : 
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(a)  Make  a  one-step  forecast 

x(t+l It)  =  ax(tit) . 


[21] 


(b)   Calculate  its  variance 


P(t+l|t)  =  a^F(tit)  +  a    . 


[22] 


(c)   Predict  the  next  observation  by 

y(t+l't)  =  x(t+l!t). 
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(d)   When  the  next  observation  is  available,  calculate  the  residual  or   "inno- 
vation" 


I(t+1)  =  y(t+l)  -  y(t+l!t) 


[24 


(e)   The  variance  of  the  innovation  is 

V(t+1)  =  P(t+1 !t)  +  R 


[25] 


;f)   If  the  errors  are  Gaussian,  the  contribution  to 
this  step  of  the  recursion  is 

ln(V(t+l)  4-  [I(t  +  1)]^  /  V(t+1). 


-2   ln( likelihood)   from 


[26] 


(g)   The  state  is  updated  by 

x(t+l!t+l)  =  x(t+l!t)  +  P(t+1 it)I(t+l)  /  V(t+1). 


[27; 


This  equation  can  also  be  written  as  a  weighted  average  of  the  old  esti- 
mate and  the  new  observation,  with  the  weighting  inversely  proportional 
to  the  variances,  giving 


x(t+l:t+l) 


x(t^lit)  ^  v(t+l) 
P(t+llt)      R 


/ 


1 


P(t+llt)    R 


[28] 


(h)   Finally  the  variance  is  updated  by 


P(t+lit+l)   =  P(t+llt)   -   [P{t-l!t)]    /  V(t+1) 


[29] 
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2 
The  -2  ln( likelihood)  function  is  calculated  for  assumed  values  of  a.  a   .     and 

R  by  summing  (26)  over  all  the  data  points: 

-2  In(likelihood)  =  S  {ln[V(t+l)]  +  [I(t+1)]^  /  V(t+1)}         [30] 

A  nonlinear  optimization  routine  then  varies  the  values  of  the  parameters 
until  the  minimum  value  for  -2  In(likelihood)  is  found.  When  data  points  are 
missing,  only  :._..ps  (a)  and  (b).  which  make  a  prediction  and  estimate  its 
variance,  are  calculated  and  then  the  algorithm  returns  to  step  (a).  Thus, 
the  above  procedure  gives  the  exact  likelihood  for  the  available  data  even 
though  observations  are  missing. 

Most  times,  the  variance  a     is  not  known:  it  can  be  removed  from  the 

2 
estimation  problem  by  differentiating  (30)  with  respect  to  a     and  setting  the 

result  equal  to  zero  (Jones . 1980) . 

a^  =  i  E  1(1)^  /  V(i)  [31] 

"  i  =  l 

Substituting  this  result  back  into  equation  (30)   gives   the  function  to  be 

minimized. 

n  n         2 

1  =  £   In  V.  +  n  In  Z   [I(t+1)]  /  V.  [32] 

1=1     ^  i=l 

Akaike's  Information  Criterium  (AIC)  is  then  computed  by 

AIC  =  -2  ln( likelihood)  +  2  [33] 

Generalizing  to  a  higher-order  model  necessitates  the  use  of  two  equa- 
tions when  using  the  Kalman  state  space  model.  The  general  case  uses  matrices 
instead  of  scalars.  The  first  equation  is  the  state  equation,  expressed  in 
matrix  form  as 


x(t+lit^l) 
x(t+2it+l) 


x(t-^mi  t+1 ) 


0    1    0 

.  .  .      0 

x(t.it) 

1 

0   0    1 

.  .  .      0 

X 

x(t+l: t) 

r 

^2 

m 

•    «2    «1 

.x(t+ni-l  It). 

.^m. 

The  state  equation  can  be  written  as 


't-il 
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[34] 


X(t  +  1)  =  (D{t)X(t)  +  U(t) 


[35] 


where  X(t)  is  the  state  vector,  <I»(t)  is  the  state  transition  matrix,  and  U(t) 
is  a  vector  of  the  random  plant  noise  with  a  covariance  matrix  Q(t).  The  ge 
term  (see  equation  (10))  corresponds  to  the  U(t)  term.  The  second  equation  in 
the  state-space  representation  is  the  observational  equation 


Y(t)  =  M(t)X(t)  +  v(t) 


[36] 


where  M(t)  is  the  matrix  showing  which  linear  combinations  of  the  state  vector 
are  observed.   The  observation  vector  Y(t)  does  not  have  to  be  the  same  length 
as  the  X(t)  vector.   The  random  observational  error  vector  v(t)   has  the 
covariance  matrix  R(t).   The  algorithm  for  the  general  case  is: 
(a)  Make  a  one-step  forecast  by 

X(t  +  ll  t)  =  a>(t)X(t!  t)  .  [37] 


(b)   Calculate  its  covariance  by 

P(t^lit)  =  «(t)P(tit)[0(t)] '  +  Q(t: 


[38] 


!c)   Predict  the  next  observation  vector  by 

Y(t+li  t)  =  M(t-l)X(t-^li  t) 


[39] 


(d)   When  the  next  observation  is  available,  calculate  the  innovation  vector 

I(t+1)  =  Y(t+1)  -  Y(t-lt)  [40] 
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(e)  The  covariance  matrix  of  the  innovation  vector  is 

V(t-H)  =  M(t  +  l)P(t-(-lit)[M(t  +  l)  j  '  -H  R{t^l)  [41] 

(f)  The   contribution   to   -2   ln( likelihood)    is 

lniV(t+i)!    +    [Kt+D] '[V(t+l)]"4(t  +  l)  [42] 

where  'V(t+l)i  is  the  determinant  of  the  covariance  matrix. 

(g)  The  state  is  updated  by 

X(t  +  lit  +  l)  =X(t  +  l't)  ^-  P(t  +  li  t)[M(t  +  l)]  '  [Vlt-fDj^^it  +  J  )  .        [43] 

(h)   Finally,  the  state  covariance  matrix  is  updated  by 

P(t  +  1  t  +  1)  =  P(t  +  llt)  -  P(t  +  l|t)[M(t)]' [V(t  +  l)]~-^M(t)P(t  +  l!t).       [44] 

The  AIC  number  is  computed  by  an  expansion  of  equation  (33): 

AIC  =  -2  In(likellhood)  +  2(p  +  q)  [45] 

The  above  procedure  proceeds  satisfactorily  as  long  as  the  coefficients 
show  no  inclination  to  pass  the  boundary  as  established  in  equations  (2)  and 
(3).  If  they  wander  close  to  the  boundary  while  being  evaluated  by  the  optim- 
ization program,  the  stability  and  invertibility  of  the  model  may  be 
threatened.  Therefore,  Jones'  subroutines  ensure  that  the  model  stays  within 
its  limits  by  reparameterizing  in  terms  of  the  partial  autoregressivo  coeffi- 
cients, a  ,  so  that  the  a's  can  be  constrained  to  lie  between  -1  and  +1  by  the 
transformation 

a^  =  [1  -  exp(-u^)]  /  [1  +  exp(-u^)]  [46] 

The  a's  can  then  be  calculated  using  the  Levinson-Durbin  recursion,   for 

n=l p.   The  n's  in  equations  (45-48)  increment  by  1,  up  to  and  including 

the  value  p.  each  time  the  recursion  is  done. 
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a   =  a  [47 1 

n    n 

If  n  =  1,  the  recursion  is  complete  and  begins  again  at  equation  (45)  with  n  = 

2.   When  n  equals  2  or  more,  a  is  calculated,  where  m  =  l,...,n-l 

m 

a  =  a  -  a  a  [48] 

m    m    n  n-m 

The  a.  variables  are  then  reset  for  the  next  recursion. 

1 

a.    =   a.  i  =  l n-1  [49] 

11  '■ 

If  n  5<=  p,  the  recursion  begins  again  at  equation  (46);  otherwise  the  autore- 
gressive  (a.)    coefficients  have  been  found  when  n  =  p. 

The  moving  average  coefficients  can  similarly  be  transformed  in  the 
recursion  (equations  (50-53)  for  n  =  l,...,q. 


bj^  =  [1  -  exp(-Wjj)]  /  [1  +  exp(-w^)],  [50] 

The  newest  ^'s  are  initialized  for  the  current  process  by 

|3   =  b  [51] 

n    n 

with  the  next  value  for  b  ,  when  n  >  1,  calculated  by 

n  ^ 

\-^m^\^m  -"^l-^ n-1  [52] 

As  in  the  autoregressive  case,  |3  =  b  when  n  =  1.   The  /3  values  are  updated 

by 


3.  =  b.         i  =  1 n-1  [53] 


and  the  process  continues  at  equation  (50)  until  n  =  q. 
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When  observation  error  is  included  in  the  model,  a  transformation  can  be 
used  to  ensure  a  nonnegative  estimate  for  the  observational  error  variance. 

2  s 

R  =  s   or  R  =  e  [54] 

Now,  the  optimization  routine  can  safely  be  run  using  the  transformed  vari- 
ables 


u,  .  k  =  1 
k 


\.    k  =  1 q 


while  the  likelihood  routine  uses  the  untransformed  variables  of  the  a's  and 

the  ^'s. 

RESULTS 

The  following  models  (Table  1)  were  selected  for  each  data  file  by  the 
criteria  of  AIC  value,  number  of  parameters,  and  suitability  to  the  present 
case.  A  pattern  to  simulate  the  data  in  each  file  can  be  generated  by  using  a 
slightly  modified  version  of  equation  (1).  Since  a  zero  mean  was  used  in  that 
equation,  the  data  simulation  series,  y,  requires  the  mean  x  value  to  be  added 
to  the  X  value. 

\  =   "  «1^-1  '   «2^-2  "   •  ••  '%\-p 

-   ^iVl  "  '^2^-2  -   •  •  •   ^  Vt-q  "  \ 

y^  =  x^  +-  X  [55] 
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The  a's  and  (3's  are  given  in  the  table  below,  as  is  the  mean.  The  e's  are 
generated  from  a  gaussian  random  number  generator  with  a  standard  deviation 
given  by  the  estimated  error  deviation  in  the  table.  For  test  files,  the 
recursive  processes  can  be  started  by  generating  q  gaussian  random  noises  and 
initializing  p  x  variables  to  zero.  The  process  then  generates  numbers.  The 
first  few  hundred  values  are  not  recorded  so  as  to  avoid  the  effects  of  the 
initial  zeros.  Any  subsequent  values  can  be  collected  and  used  to  simulate 
the  data  in  the  original  file.  A  FORTRAN  program  to  regenerate  the  data 
series  is  given  Appendix  D. 
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DISCUSSION 

Spectral  density  was  considered  as  an  analysis  technique  in  the  hope  that 
the  data  would  show  distinct  variances  at  certain  frequencies.  One  author, 
Wendenborn(1966) ,  showed  the  recovery  of  both  the  original  amplitude  and  fre- 
quency of  three  sine  functions  which  made  up  a  test  series.  The  procedure 
worked  for  that  series,  but  did  not  give  consistent  results  for  other  series. 
An  attempt  at  relating  the  amplitude  of  known  sine  curves  to  their  peaks  on 
the  spectral  density  graph  as  obtained  by  the  DQSE  method  by  regression  as  a 
function  of  original  amplitude,  length  of  the  time  series,  and  frequency  of 
the  original  curve  was  not  successful.  In  light  of  the  unpredictable  effects 
of  the  beginning  and  end  of  the  fine  data,  the  attempt  was  abandoned.  An 
attempt  to  follow  the  procedure  used  by  Cryer  and  Nawrocki  (1976)  was  impossi- 
ble because  the  original  tractor  was  no  longer  available. 

The  models  were  obtained  by  ranking  40  ARMA(p,q)  models.  The  values  of  p 
ranged  from  0  to  6  and  the  values  of  q  from  0  to  5,  with  the  sum  of  p  and  q 
less  than  or  equal  to  10,  giving  models  from  ARMA(0,1)  to  ARMA(5.5).  and 
ARMA(6.0)  to  ARiMA{6,4).  The  parameters  for  the  models  were  selected  by  a  com- 
puter program  according  to  the  the  lowest  -2  ln( likelihood) .  The  values  of 
AIC  for  each  model  were  calculated  from  the  likelihood  value  and  the  number  of 
parameters  and  the  models  were  ranked  by  ascending  order  of  AIC.  For  each 
data  file.  the  models  with  the  lowest  AIC  numbers  were  considered  in  the 
selection  of  the  best  model  of  that  data.  To  test  each  model.  a  series  was 
generated  by  setting  the  first  p  variables  to  zero  and  generating  q  random 
noises.  Values  for  x  were  then  recursively  generated  to  check  the  suitability 
of  these  models.   The  first  300  values  of  the  series  were  discarded  so  as  to 
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eliminate  the  effect  of  the  initial  zeroes.   Subsequent  values  were  then  exam- 
ined for  similarity  to  the  existing  parts  of  the  data  file. 

For  some  of  the  models  that  were  considered,  but  not  chosen,  the  gen- 
erated series  differed  from  the  original,  likely  because  the  solution  calcu- 
lated by  the  optimization  routine  was  approximate  due  to  nonlinearity  of  the 
likelihood  function  near  the  optimum  point.  These  approximate  solutions  often 
were  rated  among  the  best  by  AIC  number  and  usually  gave  reasonable  results  so 
these  models  were  considered  as  possibilities.  In  a  few  cases,  the  regen- 
erated series  oscillated  through  a  much  wider  power  range  than  was  possible  or 
was  unreasonable,  with  ridiculously  low  power  levels,  for  the  tractor. 
Another  reason  for  rejection  was  the  occurrence  of  many  values  over  the  power 
rating  of  the  tractor.  This  happened  only  when  the  mean  power  level  was  high. 
While  some  of  these  values  would  be  expected  to  occur  from  the  statistical 
standpoint  from  the  properties  of  gaussian  distribution,  for  practical  imple- 
mentation, they  would  be  undesirable.  The  third  reason  for  rejecting  a  model 
was  observation  of  repetitive  sharp  changes  in  load  level.  Although  some  of 
these  sharp  increments  were  observed  in  the  data,  for  some  numbers  generated 
by  the  models,  the  changes  were  large  and  frequent;  at  times,  changes  of  :t20 
kW  were  observed  to  occur  sequentially.  Despite  instances  where  the  parame- 
ters did  not  adequately  fit  the  data,  no  cases  were  observed  where  the 
sequence  diverged  radically  from  the  mean  and  never  changed  sign  again;  all 
the  tested  models  oscillated  about  a  mean  of  approximately  zero. 

The  models  that  remained  under  consideration  after  checking  the  regen- 
erated series  were  then  rated  by  the  closeness  of  the  AIC  numbers  and  the 
estimated  error  variance.  If  several  models  displayed  similar  AIC  and  EEV 
values,   the  simplest  model  was  chosen.   The  definition  of  "close"  was  changed 
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from  the  guideline  used  by  Jones.  The  reason  for  this  was  the  greater  value 
of  the  likelihood  from  the  data  files  used  in  this  analysis  as  compared  to 
that  of  Jones'.  If  the  models  appeared  to  give  the  same  type  of  values,  with 
widely-differing  AIC  values,  the  model  with  the  lowest  AIC  value  was  chosen. 
Some  data  files  (9a  and  10b2)  give  highly  cyclical  patterns  with  the  "best" 
choice  according  to  AIC  and  the  other  considerations.  While  the  individual 
increments  are  not  excessive,  the  series  as  a  whole  appears  to  be  more  cycli- 
cal than  the  available  data.  Thus,  alternate  models  are  offered  for  con- 
sideration. 

The  complexity  of  the  models  is  perhaps  related  to  the  amount  of  missing 
data.  In  a  test  run  with  data  generated  from  a  known  AR(3)  process,  the  AIC 
number  correctly  predicted  the  model  for  the  cases  which  used  a  complete  data 
set  and  a  pattern  of  15  values  followed  by  1  missing.  The  AR(3)  model  was  not 
selected  as  the  top  model  for  the  case  of  15  present  followed  by  5  missing  or 
in  the  case  of  15  present  followed  by  60  missing,  as  in  the  data  files.  In 
the  last  two  cases,  the  parameters  were  correct  for  the  AR(3)  model,  but  the 
AIC  number  showed  that  this  model  was  not  the  most  probable.  The  "top"  models 
were  the  more  complex  models. 

As  with  all  modeling  procedures,  roundoff  error  in  the  model-fitting  rou- 
tine is  a  concern.  Such  a  problem  was  encountered  in  6  out  of  960  models. 
According  to  Jones,  this  could  be  remedied  by  inserting  a  small  value  for 
observational  error  into  the  program,  which  has  the  option  to  consider  values 
read  with  error.  Roundoff  error  was  also  noted  when  using  a  file  of  small 
test  values  with  a  pattern  of  15  present,  60  missing.  When  the  complete  file 
was  used,  the  problem  did  not  occur.  The  problem  was  alleviated  by  multiply- 
ing the  file  by  2;  the  solution  was  adequately  obtained. 
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A  small  section  of  data  and  the  regenerated  series  is  shown  in  Appendix  C 
for  each  data  sample.  On  some  of  the  graphs  of  actual  data,  a  large  power  dip 
can  be  noted.  This  occurs  when  the  tractor  makes  a  turn  since  the  speed  is 
constant  and  the  rack  signal,  indicating  torque,  drops  sharply.  On  the 
graphs,  two  other  items  of  note  are  that  the  data  displayed  is  only  a  small 
portion  of  that  available  for  analysis  and  on  the  regenerated  series,  if  the 
power  value  exceeds  115  kW,  as  is  possible  according  to  a  gaussian  distribu- 
tion, the  graphed  value  was  115  kW,  as  it  would  be  limited  when  actually  used 
on  the  dynamometer. 

The  resulting  models  are  unique  for  each  field  data  file,  but  some  gen- 
eral characteristics  can  be  noted.  The  difference  between  models  reflects  the 
variability  between  agricultural  fields  and  operating  practices.  The  drills 
show  a  low  mean.  Both  la  and  lb  show  only  small  deviations  from  the  general 
trend,  which  is  relatively  flat.  The  2a  and  2b  series  show  a  rougher  pattern, 
particularly  2a.  While  some  these  changes  may  appear  to  be  excessive,  the 
actual  does  show  some  large  steps.  Sample  la  is  chosed  to  represent  the 
drills  since  the  actual  and  regenerated  series  are  quite  similar  and  resemble 
the  rest  of  the  samples. 

The  series  for  a  chisel  (location  4b)  is  quite  rough.  This  is  expected 
since  this  tool  is  a  primary  tillage  tool  with  a  fracturing  mode  of  soil  frac- 
ture as  opposed  to  cutting  by  a  disk.  The  mean  level  is  relatively  high  (79.7 
kW)  and  the  series  show  power  fluctuations  reaching  almost  to  the  maximum 
power  level. 

The  drag/harrow  combination  (locations  7a2  and  7b2)  showed  a  smooth  pat- 
tern.  The  simulated  series  for  7b2  does  not  show  the  sharp  peaks  that  some  of 
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the  data  shows.  The  drag,  as  recorded  by  Blumanhourst ,  was  an  old  truck 
frame.  This  combination  was  not  hooked  soiidly  to  the  tractor  and  so  would 
tend  to  jerk  at  the  corners.  The  smoothness  of  the  series  is  expected  since 
this  tool  combination  was  being  used  to  level  and  smooth  the  seedbed  under  a 
center  pivot  irrigation  system.  The  operation  can  be  represented  by  location 
7a2. 

The  disk/springtooth  combination  (locations  7al  and  7bl )  shows  a  pattern 
of  small  oscillations  about  the  long-term  trend.  The  data  for  series  7  shows 
the  effect  of  turns  with  this  operator.  The  combination  of  tools  does  not 
appear  to  be  distinct  from  the  disk  alone,  given  the  variability  of  disk  sam- 
ples. The  series  are  quite  similar,  but  since  7al  avoids  the  apparent  cycli- 
cal tendency  of  7bl,  an  ordinary  field  without  definite  slopes  can  be 
represented  by  the  7al  series. 

The  field  cultivator  (locations  9b2  ,  10b2,  and  10b3)  provide  the  oppor- 
tunity to  observe  the  differences  between  operators  and  fields.  Location  9b2 
had  an  actual  mean  of  62.9  kW,  and  10b2  is  89.1  kW.  while  at  10b3  the  mean  was 
96.3  kW.  As  a  result,  10b3  shows  a  series  with  some  oscillations  near  the 
maximum  level  while  9b2  never  approaches  these  levels.  For  9b2,  the  simulated 
series  varies  from  the  mean  with  no  evident  pattern  while  10b2  and  10b3  have 
some  cyclical  tendencies.  At  location  10,  the  apparent  small  field  size,  with 
the  resulting  frequent  drops  in  power  levels  at  the  turns,  may  have  contri- 
buted to  a  cyclical  tendency  and  large  step  size  between  the  proposed  patterns 
for  10b2  and  lObS  while  the  data  shows  a  finer  variation.  Thus,  9b2  is  the 
best  choice  to  represent  a  field  cultivator.  At  10b2,  the  alternate  choice 
agrees  more  with  the  series  for  3b.  which  is  NH3  application,  in  that  the 
variations  from  one  time  increment  to  the  next  are  smaller.   The  ground 
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surface  would  be  expected  to  be  fairly  even  when  applying  anhydrous  ammonia 
since  the  soil  surface  needs  to  be  fine  enough  so  the  ammonia  does  not  escape 
from  the  soil . 

The  disk  tool  has  the  largest  number  of  samples  and  also  the  largest 
variation  between  the  eleven  samples.  Four  sizes  of  disks  are  represented. 
5.5m,  6.1m,  6.7m.  and  9.1m.  The  mean  power  levels  range  from  39.8  kW  to  88.3 
kW.  Most  samples  and  regenerated  series  show  a  moderate  variation  from  one 
Increment  to  the  next.  The  long  term  trend  for  some  is  cyclical  (6a  or  9a, 
alt.).  For  others,  the  trend  is  flat  (8a  or  8b).  Size  of  tool  does  not 
appear  to  be  a  determining  factor  for  the  variation;  the  EED  for  9b  (5.5m 
disk)  is  .923  kW,  for  6b  (9.1m  disk)  it  is  1.086  kW.  and  for  8b  (6.7m  disk) 
the  deviation  is  .832  kW.  There  is  also  variation  between  samples  for  the 
same  location.  For  location  8,  8a  has  1.826  kW  EED,  8b  has  .832  kW,  and  8blo 
has  1.790  kW,  yet  the  means  for  8a  and  8b  are  close.  Given  the  variability 
between  each  series,  a  "typical"  disk  pattern  is  difficult  to  select.  In  the 
interests  of  providing  a  pattern  intermediate  between  the  smooth  drag/harrow 
and  the  chisel,  which  most  of  the  samples  appear  to  be,  8a  is  selected  because 
it  contains  both  small  incremental  changes  and  sharp  adjustments  in  power 
usage. 

CONCLUSIONS 


1.  The  results  of  this  study  are  a  model  of  each  of  the  data  samples.  The 
models  were  obtained  using  the  ARMA(p,q)  approach  in  which  the  value  at 
time  t  is  dependent  on  p  past  values,  q  past  random  noises,  and  a  random 
noise  at  time  t. 

2.  Generalizations  about  the  models  are  possible.  Drills  and  the 
drag/harrow  combinations  have  a  relatively  smooth  pattern.  The  chisel 
showed  a  rough  pattern.  The  field  cultivator  has  a  roughness  less  than 
the  chisel,  but  more  than  the  disk.   The  disk  samples  showed  a  variety  of 
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patterns,  but  the  power  loading  pattern  for  each  model  was  generally 
between  the  drill  and  chisel  in  terms  of  roughness.  The  disk/springtooth 
combination  was  not  noticeably  different  from  the  disk.  The  variation 
when  NH3  was  applied  was  apt  to  be  small  between  successive  time  incre- 
ments . 


SUMMARY 

The  autoregressive  moving  average  approach  was  used  to  determine  a  model 
for  each  data  sample  which  could  recursively  regenerate  the  represented  series 
of  power  loading  variations  for  a  tractor.  The  method  accommodated  evenly- 
spaced  data  with  missing  values  by  adjusting  coefficients  for  the  past  values 
through  a  nonlinear  optimization  routine  until  the  minimum  value  of  -2 
ln( likelihood)  was  reached.  Spectral  density  was  also  considered,  but  was 
rejected  due  to  the  difficulty  of  definitely  correlating  the  amplitude  of 
known  sine  curves  to  the  values  on  the  power  spectral  density  graph. 

The  models  of  power  loading  variation  obtained  for  the  different  tillage 
tools  show  that  the  drill  and  drag/harrow  have  fairly  small  variations  between 
time  increments.  The  disk  samples  generally  show  a  moderately  rougher  pat- 
tern, with  no  noticeable  difference  when  a  sprlngtooth  is  attached.  The 
chisel  appears  to  be  the  roughest  series  and  the  field  cultivator  shows 
slightly  smaller  deviations.  When  NH3  is  applied,  the  increments  from  step  to 
step  are  small  although  the  general  trend  varies. 

SUGGESTIONS  FOR  FURTHER  RESEARCH 

Anyone  wishing  to  do  work  on  the  time  series  of  tractor  power  variations 
with  various  tillage  implements,  or  any  other  time  series,  would  find  it  bene- 
ficial to  begin  with  a  series  with  no  missing  values  at  evenly-spaced 
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intervals.   This  would  allow  the  use  of  conventional  theory  for  which  means  of 
analysis  have  been  well-established. 

The  time  series  analysis  in  this  study  was  done  on  field  data  without 
quantifying  the  effects  of  different  variables  on  the  power  level.  A  careful 
analysis  could  lead  to  an  even  better  representation  of  power  usage  during 
tillage.   Possible  areas  to  investigate  include  the  following: 

1.  The  effect  of  different  soil  types  on  the  power  usage  pattern.  Soils  of 
varying  composition  could  possibly  yield  distinct  patterns. 

2.  The  differences  between  the  sequences  derived  from  specific  tillage 
implements . 

3.  The  impact  on  the  regenerated  series  of  the  actual  power  drops  due  to 
turns.  When  the  tractor  makes  a  turn,  the  power  level  drops  sharply  and 
then  rises  sharply.  The  regenerated  series  was  not  able  to  show  these 
sharp  drops . 

4.  Simulation  of  hills  with  the  possibility  of  adding  or  removing  slopes 
when  desired. 
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APPENDIX  A  -  CONTINUOUS  AND  DISCRETE  SPECTRA 

Continuous  Spectra 

The  spectral  density  method  of  analyzing  time  series  relates  the  autoco- 
variance  to  the  frequency.  A  plot  of  power  spectral  density  against  frequency 
shows  whether  a  series  is  white  noise,  with  a  constant  variance  over  all  fre- 
quencies, or  has  a  periodic  tendency  within  a  certain  frequency  range.  The 
main  sources  for  this  presentation  are  Blackman  and  Tukey  (1958)  and  Otnes  and 
Enochson  (1972).  The  theory  can  also  be  found  in  other  sources  dealing  with 
time  series  analysis. 

To  begin  with,  let  the  time  series  x(t)  be  a  simple  sinusoid. 

X(t)  =  cos  (tit  [1] 


For  simplicity,  let  w  be  a  nonnegative  number.  This  series  goes  through  one 
complete  cycle,  a  period,  in  time  T.  The  frequency,  f,  is  the  reciprocal  of  T 
and  is  measured  in  cycles/sec  (hz).  So  the  units  cancel,  w  is  defined  as  the 
angular  frequency  (rad/sec).   By  definition,  2k   radians  equal  a  cycle  so 

w   =  27rf  [2] 


Adding  two  more  parameters,  A  and  0,  called  amplitude  and  phase,  respectively, 
gives  a  more  general  equation 

x(t)  =  A  cos(wt  +  (p)  [3] 


A  is  the  maximum  amplitude  from  the  mean  value,  in  this  case  zero,  and  <p  is 
the  displacement  of  the  sinusoid  from  the  given  time  origin  (the  phase  shift). 
For  simplicity,  A  is  a  non-negative  number.  Any  number  of  patterns  can  be 
created  by  changing  the  parameters  in  equation  (3)  and  adding  one  or  more 
sinusoids  together.  Most  natural  processes  would  be  expected  to  contain  more 
than  one  frequency.   This  new  type  of  series  is  given  by 

n 

x(t)  =   I  A.cos(cw.t  +  0.).  -00  <  t  <  0°  [4] 

.,111  ■• 

1  =  1 

The  classic  approach  to  the  power  spectral  density  as  outlined  by  Blackman  and 
Tukey  (1958,  pp. 84-87)  deals  first  with  autocovariance .  Let  x(t)  be  some  real 
function  of  time.  If  the  mean  of  the  function  is  0  as  time  approaches  time 
limit  T 

T/2 

rtl  ?  ^  ><(t)  dt  =  0  [5] 

-T/2 
then  the  autocovariance  function,  R,  of  the  process  is  defined  as: 
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T/2 
•^(^^  =  1^1  T       i       x(t).x(t^r)  dt  [6J 

-T/2 

The  variable  t  is  a  difference  in  time  known  as  the  time  lag  and  is  unre- 
lated to  (t>.  Where  R(t)  is  defined  for  a  finite  piece  of  a  function,  it  is 
called  the  apparent  autocovariance  function.  When  t  is  0,  a  special  case  of 
the  autocovariance  function  results,  the  variance  function: 

T/2 
R(0)  =111^       I   [x(t)]^  dt  [7J 

-T/2 

This  definition  of  variance  is  the  same  as  the  usual   statistical  definition 
since  the  function  x(t)  has  the  mean  defined  to  be  zero. 

The  spectral  density  function  is  the  Fourier  transform  of  the  autocovari- 
ance density. 

P(f)  =  I  R(r).e"^"''^dr  [8] 


For  a  derivation  of  this  as  well  as  an  explanation  of  the  use  of  the  word 
"power"  to  describe  the  spectrum,  see  Blackman  and  Tukey  (pp. 85-86).  K(f)  is 
the  inverse  Fourier  transform  of  P(r): 

00 

R(r)  =  I  P(f)  e^^'^df  [9] 


Even  and  odd  functions  have  properties   that  can  simplify  an  equation 
(Boyce   and  DiPrima,  1977).   P(f)  and  R(r)  are  even  functions  since  P(f)  =  P(- 
f )  and  R(r)  =  R(-r) . 
Letting  g(T)  be  an  even  function  as  shown  above, 

00  09  0 

;  g(z)dv   =   J  g(r)dr  +  J  g(r)dT 
substituting  r  =  -s  in  the  rightmost  term, 

00  0 

=1  g{r)d(r)  -  ;  g(s)  ds 

0  oo 

The  -  sign  comes  from  the  derivative  of  s  [g(s)  =  g(-s)].  Inverting  limits  on 
the  rightmost  term  adds  another  -  sign,  making  the  overall  sign  of  the  term 
positive . 

00  00 

=  I  g(r)dr  +  J  g(s)ds 
0  0 
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=  2  ;  g(T)dr 
0 

The  area  in  both  integrals  is  the  same  and  additive. 

When  h{r)  is  an  odd  function,  h(-r)  =  -h(r)   and  the  integration  from 
-00  to  +■»  gives  a  result  of  zero. 

00  00  0 

J  h(r)dr  =  J  h(T)dr  +  J  h(T)dr 

-oo  0  -<" 

oo  0 

=  J  h(r)dr  +  j  -h(s)(-ds) 
0  " 

00  CO 

=  ih(r)dr  -  j  h(s)ds  =  0 
0         0 

The  area  in  both  integrals  is  the  same,  leading  to  an  answer  of  0. 
Equation  (8)  can  be  rewritten  as 

00 

P(f)  =  J  R(t)[cos(ut  -  i  sin  u)z](iz  [10] 

The  function  cos  wz   is  even  since  cos(«>r)  =  cos(-u)z).        Multiplication  of  an 
even  function  by  another  even  function  yields  an  even  function: 
Let  g(r)  =  R{r)cos(a)T) 

g(-r)  =  R(-r) -cosl-wr) 

=  R(r)  •  cos{<ur ) 

=  g(r) 

Multiplication  of  an  even  function  by  an  odd  function  yields  an  odd  function. 
The  sin  (wr)  is  an  odd  function  since  sin(a)T)  =  -sin(a)r). 
Let  h(r)  =  R(T)sin  cur 

h(-r)  =R(-r)sin(-a)T) 
=  R(r ) [-sin  wz] 
=   -R(r)sin  wz 
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-h(T) 


Since  the  integral  of  an  odd  function  h(r)  is  zero,  equation  (10)  can  be  sim- 
plified to  the  cosine  transform  of  R(r) 


P(f)  =  I  R(t)cos  o)z   dr 


and 


R(z)    =  I  P(f)cos  2KfT   df 


[11] 


(2Kf     =  (O) 

or  even  more  simply  as  one-sided  cosine  transforms  using  the  property  of  an 
even  function: 


P(f)  =  2  I  R{r)cos  uiz   dr 
0 


[12] 


R(r)  =  2  J  P(f)cos  2;rfr  df 
0 


[13] 


For  further  discussion  of  the  power  spectral  density,  two  more  mathematical 
concepts  will  be  helpful.  One  of  these  is  known  as  the  Dirac  delta  function. 
The  delta  function,  8(t-t  )  is  introduced  by  formally  identifying  f(t  -  ty)<lt 
with  dh(t  -  t  )  (BlacSman  and  Tukey,  1958, p.  69)  where  h(t  -  t^)  is  Heavi- 
sides'  step  function 


h(t 
h(t 


t„)  =  0, 
^o^  -  ' 


t  <  t 
c 
t  <  t 


So  if  the  time  function  is  the  delta  function 


G(t)  =  6(t-t  ) 


the  equivalent  frequency  function  is 


S(f) 


-iwt 
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The  delta  function  is  arbitrarily  tall  and  narrow,  yet  has  a  unit  area.  Otnes 
and  Enochson  (1972,  p. 14)  also  give  two  functions  derived  from  the  delta  func- 
tion. 


■1 


5(f-f  )  -  5(f-f„, 
o         o 


cos  27:f  t 
o 


sin  2;rf  t 
o 


[15] 


[16] 
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Convolution,  another  important  property,  is  the  operational  transform 
between  the  frequency  and  time  domain  that  works  as  follows  (Blackman  and 
Tukey,1958,  p. 72):  If  G(t)  =  G  (t)  '  G  (t),  and  S(f)  is  the  frequency  counter- 
part of  G(t),  then  the  Fourier  transform  of  G{t)  in  terms  of  G  (t)  and  G  (t) 
is 


S(f)  =  I  G^{t).G2(t).e 


-i(ot 


dt 


[17] 


=  ; 


G^(t) 


J  S^in.e^'^^^dl 


•  e    dt , 


[18] 


where  S  (f)  is  the  frequency  domain  equivalent  of  G  (t) 


=   J 


;  G^(t).e-i2.(f-e)t^^ 


S2(€)d§ 


[19] 


=  J  S^(f-e)S2(§)d^ 


[20] 


S  (f)  and  S  (f)  are  interchangeable  in  this  operation,  also  written  as 


S(f)  =  Sj(f)  *  S^ff) 


[21] 


The  "*"  indicates  the  convolution  operation.  For  the  completeness  of  the 
transform  operation,  the  process  of  going  from  multiplication  in  the  frequency 
domain  to  convolution  in  the  time  domain  is,  given  S(f)  =  S  (f)  '  S  (f). 


G(t)  =  J  G^(r-/5)  .  (^^KX)AX   =  G^(t)  *  G^tt) 


[22] 


For  the  continuous  domain,  the  limits  of  +«>  to  -«>  work  well,  but  In  the 
real  world  of  experiments,  the  data  samples  have  a  beginning  and  an  end  as 
well  as  a  finite  length.  To  deal  with  these  difficulties,  the  idea  of  a  "box- 
car function"  is  introduced,  designated  as  u  (r). 


u.^  (  r )  =  1    I  r  I  ^  T 


u^{r)  =  0 


>  T 


Its  Fourier  transform  is 


U^(f)  =  I  u^(t)e  ^"'^dt 

—  00 


T 

J  e 
-T 


-io)t  _  2  sin  oT 


[23] 


The  infinite  limited  transform  integral  for  the  power  spectrum  can  now  be 
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written  in  terms  of  the  boxcar  function  (Otnes  and  Enochson,  1972).  The  prime 
sign  shows  the  variable  is  computed  for  the  finite  length  rather  than  the 
infinite. 

00 

P'(f)  =  I  u^(T)  .  R{r)  -  e""''"^dr  [24] 

—  00 

As  shown  above  in  the  convolution  theorem,  multiplication  in  the  time  domain 
corresponds  to  convolution  in  the  frequency  domain.  Thus,  equation  (24)  can 
be  written  as 

P'^m  =  P^(f)  *  U^(f)  [25] 

00 

=  J  S(fJ.U,j,(f-f^)df^  [26] 

—  00 

Substituting  for  the  U,j,(f-f^)  term  (see  Equation  (23))  yields 

sin  [2;rT(f-f  )] 
P'x(f)  =  ^  ^x^^o^  2     2.(f-f  )°   '^^o  t^^^ 

-00  '      O 

As  can  be  seen,  the  frequency  components  are  not  dependent  on  one  another  and 
hence  their  contributions  to  the  power  spectrum  are  additive;  each  source  con- 
tributes separately. 

For  an  example  (Otnes  and  Enochson,  1972,  p. 202),  let  the  K  (r)   function 
be  the  covariance  of  a  cyclical  batch  of  data:  ^ 

2 
R  (r)  =  ^  cos  2;tf  z  [28] 

Then,  R^(f)T)  corresponds  to  P^H)    in  the  frequency  domain  (see  equation  (15)) 

P^(f)  =  ~-    [<5(f-f  )  +  5(f  +  f  )],  [29] 

X       4        o         o 

which  are  delta  functions  centered  at  -f  and  -^f  .   When  this  value  for  P  (f) 
is  substituted  into  Equation  (27)  the  result  for°the  truncated  P'(f)  is   ^ 

.2  sin[27tT(f-f^)]    2  sin[2;rT(f^f  )] 

pi  I  f  \     _    2 O       A 0 

X*  '    2     2;t(f-f  )       2.    2;r(f-f  )  ^30] 

o  o 

For  positive  frequencies  where 

0  «  i  «  f^_ 
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P',^(f) 


/I   sin[2;rT(f-f  )1 
A o  ' 

2     271  (f-f  ) 

0 


[31] 


The  maximum  P'  (f)  is  reached  when  the  denominator  approaches  zero 
(f-»f  ).  Since  an'  indeterminant  form  is  reached  (0/0)  when  the  limit  is 
applied,  the  maximum  P'  (f)  is  determined  by  L'Hospital's  rule  to  be 


P'  (f) 
X   ' 


lim 
f-f 


9  --  [sin  27rT(f-f  )J 
A  dt  o 


df  ^^^'-'o^ 


[32] 


,  .  ,2  2;tT  cos  27rT(f-f  )     2^ 

lira  A_  o   _  A  T 

f-f  2         27t        '2 

0 


[33; 


As  T  gets   larger,   the  peak  becomes  higher  and 


narrower,   until   T 
Since  T  is  fin- 


approaches  infinity,  when  P'  (f)  approaches  a  delta  function, 
ite,  what  would  have  been  a  very  sharp  peak  with  the  power  concentrated  within 
a  very  narrow  frequency  range  now  has  become  a  shorter,  broader  curve  with  the 
power  spread  out  due  to  the  (sin  x)/x  term  as  shown  in  Figure  (Al).  This 
spread  of  power  is  termed  leakage. 


PSD 


Figure  Al .   Boxcar  function  window 

The  first  zero  crossing  to  the  right  will  occur  at  1/4  of  the  cycle. 

27rT(f   -  f 

f 
r    o   4J 

f 


o 

4T 


2 


r    o 
For  the  left  side,  the  first  zero  crossing  occurs  at  f 


of  the  curve  is  thus  due  only  to  the  length  of  data  sample: 


4T' 


^  =  ^o  "  H 


[f 


4T^    2T 


The  spread 

[34] 
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The  areas  to  the  side  of  the  main  peak,  referred  to  as  side  lobes,  may  be 
of  a  significant  size  and  cause  the  main  peaks  to  be  distorted,  particularly 
if  frequencies  in  the  data  are  located  close  together  and  the  side  lobes  hap- 
pen to  be  in  phase  so  that  they  are  additive.  The  distinct  beginning  and  end 
of  the  data  have  the  effect  of  inserting  a  (sin  x)/x  term.  The  spreading  can- 
not be  eliminated,  but  different  windows,  as  the  the  type  of  function  (sin 
x)/x  is  referred  to,  may  be  used  to  increase  the  height  of  the  main  peak  or 
increase  the  spread  of  the  main  peak.  These  windows  can  operate  directly  on 
the  data  in  time  domain  or  on  the  transformed  values  in  the  frequency  domain. 
A  data  window  operates  in  the  time  domain  for  a  given  interval  with  the  effect 
of  multiplying  data  or  signals  that  are  defined  for  a  longer  period  and  which 
will  then  be  subjected  to  further  processing. 

Windows  may  seem  to  be  a  complicating  addition,  but  they  are  unavoidable. 
By  default,  the  window  present  is  the  basic  boxcar,  which  has  some  undesirable 
characteristics  that  can  be  reduced  by  other  windows.  These  windows  tend  to 
smooth  the  data  for  improved  quality  in  the  frequency  analysis.  The 
categories  of  lag  and  spectral  windows  appear  most  often  in  literature.  The 
lag  windows  are  a  function  of  lag,  a  time  difference  between  two  events  con- 
sidered together.  These  windows  exist  for  a  certain  interval  and  then  vanish. 
Their  counterparts  in  the  frequency  domain  are  the  spectral  windows.  Blackman 
and  Tukey  ( 1958 , pp .95-99)  list  five  possible  lag-spectral  window  pairs.  The 
boxcar  window,  inherent  in  finite  data  was  discussed  above.  Two  windows  that 
were  considered  in  the  data  analysis  were  the  banning  and  hamming  window. 

The  banning  window  is  defined  by 

=0  r  <  -T 


Uj     (T)  -  -  (1  +  cos  — ) 

m  m 


-T   ^  r  ^  T 
m       m 


[35] 


=  0 


r  >  T 


T   is  the  length  of  autocorrelation  used. 
The  Fourier  transform  of  this  is 


u^  (f)  =  K  <f)  ^  f 
m  m 


sin    [T    {2Kf   -  :^)\  sin[T    (2Kf    +■  i^) 
ml  ml 
m  m 

+  


TJ2Kf    -   ^) 
m  T 

m 


TJ2k{   +   - 
m  T 


[36] 


i  I' 
2    ^T 


f)  -  i^T  <f  -  ^'  -  i^T  (f  ^  if' 

m  m  m  m  m 


[37 


Thus  the  banning  window  turns  out  to  be  the  summation  of  three  of  the  boxcar 
functions.  Since  these  lobes  are  spaced  1/(2T  )  Hz  apart,  tbe  banning  window 
has  the  effect  of  averaging  three  adjacent  lobes  and  giving  the  center  one 
twice  as  much  weight  as  the  end  pair,  reducing  the  center  lobe  to  1/2  of  its 
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former  height,  but  also  doubling  its  width.   The  distance  between  the  first 

zero  crossings  on  either  side  of  the  main  lobe  is  now  2/T  while  the  effective 

bandwidth,  the  frequency  interval  between  half  power  poinds,  is  now  1 /T 

m 

The  second  spectral  window  is  the  hamming  window,  defined  by 
=0  z   <   -T       m 

u   (r)  =  .54  +  .46  cos  f^      -T   <  r  <  T  r-^„-, 

T  T         m        III  L>>oJ 

m  m 

=  0  r  >  T 

m 

After   Fourier   transform, 

U^    (f)    =    .54    L-^    (f)    +    .23   U^    (f    ^    i^)    ^    -23   U^    (f    -   ^)  ^39^ 

ID  m  mm  mm 

So,  the  hamming  window  is  similar  to  the  banning  window  in  that  they  are  both 
sums  of  (sin  x)/x  terms,  but  with  different  weights.  Two  differences  between 
them  can  be  noted: 

1.  The  height  of  the  maximum  side  lobe  for  the  hamming  window  is  approxi- 
mately 1/5  that  of  the  banning  window. 

2.  The  heights  of  the  side  lobes  of  the  banning  window  drop  more   rapidly 
than  do  those  of  the  hamming  window. 

Thus  one  difference  favors  the  hamming  and  the  other  favors  the  banning  win- 
dow. 

For  either  the  banning  or  hamming  window,  the  adjacent  side  lobes  with 
their  positive  and  negative  values  will  tend  to  cancel.  If  the  spectrum  has  a 
strong  power  level  at  a  distance  which  is  resonant  with  a  low  main  peak  and 
the  negative  side  lobes  tend  to  cancel  it  out,  the  low  peak  can  still  be  noted 
by  the  presence  of  the  side  lobes.  Confidence  intervals  for  spectral  densi- 
ties can  be  obtained  by  using  both  a  positive  window  such  as  Parzen's  and  a 
negative  window  (See  Wonnacott  (1961)). 

Discrete  Spectra 

The  popularity  of  digital  data  collection  indicates  that  discrete  inter- 
vals may  be  of  more  use.  Accordingly,  the  equations  are  modified  to  handle 
discrete  data  by  using  summations  rather  than  intcjgrals  and  At  instead  of  dt. 
The  following  approach,  which  has  been  described  by  Blackman  and  Tukey  (1958) 
and  Otnes  and  Enochson  (1972  p. 270),  has  been  used  to  compute  the  power  spec- 
tral density  (Walls,  et  al.,  1954)  and  (Wendenborn,  1966). 

1.  Let  X  be  a  sequence  of  N  evenly-spaced  points  having  a  zero  mean;  if  this 
is  not  the  case,  find  the  mean  and  use  the  set  of  standard  deviations. 
Wendenborn  (1966)  and  Walls,  et  al .   (1954)   use   a  running  average   to 
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eliminate   the   low  frequency  cycles   that  would  not  be  adequately 
represented  in  the  data. 


2.  The  sample  autocorrelation  R  is  calculated  for  (m+1)  values  where  m 
represents  the  maximum  "time"  lag.  Since  the  data  is  evenly  spaced,  the 
total  time  represented  is  mAt ,  so  the  number  lag  becomes  a  time  lag. 

N-r 

r  =  0,1, .. .m  [40] 


1 


N-r  .  ,   1  i+r 
1  =  1 


The  power  spectral  density  is  calculated  for  various   frequencies,   using 
trapezoidal  integration. 


r^  '    X 


m-1 


;trq 


R   +  2   I  R   cos  - — ^  -  R  cos  ;rr 
o       ,   q      mm 
q=l 


[41] 


These  P'  (f)  values  are  at  evenly-spaced  frequencies   (rad/unit  time) 
given  by 


f  = 


m  At 


r  =  0.1 


[42] 


These  raw  spectral  density  estimates  are  then  smoothed  using  a  spectral 
window.   For  example,  if  the  hamming  window  is  used: 


P  "(f )  =  .23  P' ,   ,  ,(f )  +  .54  P' (f ) 
r  (q-1)  q 


23  P'      ( f  ) 


a  ;<=  0,m 


[43] 


P-^lf)  =  .54  P'o(f)  +  .46  P',{f) 


[44] 


P"  (f) 
m 


.46  P'_^_j(f) 


.54  P'  (f) 
m 


[45] 


Step  3  could  be  done  by  alternate  methods,  such  as  the  Fast  Fourier  Transform, 
that  would  decrease  the  computer  time  needed  if  the  data  is  evenly  spaced. 
For  equally-spaced  data,  the  frequency  range  of  the  estimated  spectrum  ranges 
from  0  to  l/(2At),  the  Nyquist  folding  frequency.  Any  frequency  present  in 
the  data  that  is  higher  than  the  folding  frequency  will  be  represented,  by 
aliasing,  in  the  principal  part  of  the  spectrum.  When  0  <  f  <  l/2At,  f  is 
called  the  principal  alias  and  could  possibly  represent  its  aliases  f, 

etc . .  with 


2f  4-f,  4f  - 
N   '    N 


f.  4fj^.f. 


f   equal  to  the  Nyquist  frequency. 


2f  -f. 
N 


The  discrete  method  was  the  first  attempt  to  analyze  the  data.  Indivi- 
dual minutes  of  fine  data  were  analyzed  and  the  resulting  spectra  were  aver- 
aged to  obtain  an  "overall"  spectra.  Sequences  of  coarse  data  (minute  aver- 
ages) were  also  analyzed  for  low  frequency  patterns.  The  results  from  the 
second  analysis  were  not  practical  for  a  regenerated  data  cycle  of  10-15 
minutes  and  the  results  from  analyzing  just  15  points  at  a  time  were  incon- 
clusive so  a  longer  series  was  proposed  by  recombining  the  data.  The  litera- 
ture listed  no  method  to  recombine  data  with  unusual  spacing  so  Blackman- 
Tukey's  approach  for  evenly-spaced  discrete  data  was  abandoned. 
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APPENDIX  B  -  DIRECT  QUADRATIC  SPECTRUM  ESTIMATION 

Ways  of  handling  a  real  data  sample  for  the  discrete  evenly-spaced  case 
have  been  discussed  in  Appendix  A.  Further  work  on  time  series  has  generated 
methods  of  dealing  with  gaps  in  data  and  even  irregularly-spaced  data.  One 
method  for  dealing  with  irregular-spaced  data  is  the  Direct  Quadratic  Spectrum 
Estimation  (DQSE)  method  developed  by  Marquardt  and  Acuff.  This  appendix  sum- 
marizes the  DQSE  method  (Marquardt  and  Acuff.  1984)  which  was  used  to  deter- 
mine the  spectrum  for  the  data  on  an  irregular-space  basis. 

Several  assumptions  are  necessary  to  apply  the  DQSE  method.  It  is 
assumed  that  the  pattern  of  irregular  spacing  is  not  related  to  the  stochastic 
properties  of  the  process  y(t)  being  sampled,  where  the  independent  variable  t 
represents  either  time  or  distance  units.  The  process  y(t)  is  assumed  to  be 
at  least  weakly  stationary,  meaning  that  the  covariance  between  two  points  t. 
and  t.,  where  t.  <  t.,  depends  only  on  the  time  difference  (t.  -  t.).  The 
mean  ii^  assumed  to  be  zero.  The  data  comes  from  a  process  y(t)  in  which 
observations  are  made  at  times  t.  during  a  single,  continuous  period  of  obser- 
vation [O.Tl.   Data  is  in  the  form  (t.,Y.),  where  i=l,2,3,...n  and  0  ^  t.  ^  T. 

The  maximum  frequency  range  of  the  spectrum,  if  evenly  spaced,  is  from  0 
to  l/(2At);  any  frequency  in  the  data  which  is  higher  than  the  folding  fre- 
quency is  represented  by  its  alias.  For  irregularly  spaced  data,  the  cutoff 
point  for  folding  is  less  sharp,  particularly  so  as  the  spacing  becomes  more 
random.  The  sensitivity  to  high  frequency  cycles  below  the  Nyquist  frequency 
is  diminished,  but  some  cycles  above  the  cutoff  frequency  can  be  obtained. 
Jones  (1962)  provides  a  theoretical  base  for  the  direct  quadratic  spectrum 
estimation  using  irregular  data. 

The  DQSE  equation  takes  the  form 

n  n 
P(a,)  =  ;f7  ^  5:  W  Y(t.)Y(t  )  [1] 

i<j 

T'    is  the  effective  record  length.  Depending  on  the  data  spacing,  it  is 
equal  to  the  record  length  or  less  (see  equation  (8)). 

W.  .    is  the  weight  represented  by  the  l,j  element  in  the 
symmetric  quadratic  form  (W..  =  W..) 

The  weight  element  is  dependent  on  several  factors: 

W.  .  =  D(T:)F(t.  ,t  .)  cos  2k(i).z  [2] 

ij  1   J 


D(r)      is  the  lag  window  for  lag  r,    r  =  t.-t. 

F(t,.  ,t_.)   is  the  "data  spacing  factor",  more  precisely  defined 

by  equation  (9  or  10) 

frequency  in  cycles  per  unit  time  or  distance 


1   J 
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The  window  could  be  one  of  those  of  discussed  in  Appendix  A  or  another  suit- 
able window;  the  banning  window  was  used  by  Marquardt  and  Acuff.  The  hamming 
window  gave  similar  results  in  the  present  analysis.  So  the  window  takes  the 
form 


1  r 

D(t)  =  -  [1  +  cos  K   — ] 

m 


z   ^   r 


[3] 


=  0 


T  >  r 


[4] 


T  is  the  maximum  time  (or  distance)  lag  used 


The  data  spacing  factor  is  based  on  the  idea  that  the  weight  given  to 
each  ij  element  should  be  proportional  to  the  time  represented  by  that  ele- 
ment. The  spacing  factor  was  suggested  in  work  done  by  Jones  (1962).  In 
terms  of  each  i,j  element,  the  data  spacing  factor  is 


F(t, 


^J^  = 


t.   +t. 
1  +  1   a 


t.+t. 
1   1- 


t  .   +t  . 


t.+t.  , 


[5] 


[t, 


i(  +  ) 


'i(- 


"'ji 


-'j(-)' 


=  Vj 


The  endpoint  conditions  are  t 

the  time 

maximum  length  used  for  defining  the  sinusoid,  that  corresponds  to  the  fre 

quency.   Thus,  an  upper  bound  is  placed  on  the  weight  given  to  an  area  by  lim- 


^_j  ^  0  and  t.   ^  T.   For  large  time  increments, 
intervals  may  be  greater  than  174  of  the  cycle  length,  which  is  the 


itmg  t. 


i(  +  ) 


^i(-)- 


t.,  ^,  and  t.,  ,  where  necessary, 
J(+)       J{-)  * 


where 


t. 


i(  +  ) 


=  mm 


t.  ,  +  t. 


^'i 


40)' 


[6] 


^i(-) 


^  max 


t . 

,  1 


^i-l 


), 


:^i  -  ^^ 


The  <5  factor  is  expressed  similarly  by  substituting  j  in  place  of  i  on  the 
right  side  of  the  last  two  equations.  End  conditions  are  defined  at  t  =  0  and 
t  =  T  by  setting  t^  =  -t^ ,  and  t  =  2T-t  .  When  (o  =  0,  or  is  still  small, 
the  periods  are  very  large,  with  the  possible  result  of  large  data  gaps  being 
within  the  minimum  5,  and  having  an  undue  influence  on  the  spectrum.  The 
value  of  (1)   is  thus  constrained  to  be 


a>(see  equation  6)  =  max[(u,  Mo  J 

M 


[7J 
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The  variable,  M,  is  an  arbitrary  multiplier  (.1  is  given  as  an  example)  and 
(t)  =  (n-l)/2T)  is  the  Nyquist  frequency  associated  with  the  average  inter- 
point  interval  in  the  given  record. 

Equation  (6)  implies  a  decrease  in  the  effective  length  of  time  when  the 
data  is  irregularly  spaced  or  contains  gaps.  The  effective  length  of 
T' (w)  ^  T  is  defined: 

n 

T'  ((1))    =  Z  (5.  [8] 

i  =  l 

The  (5.  are  defined  in  equation  (6).  If  the  data  is  evenly  spaced  with  no 
gaps,  the  effective  length  simplies  to  the  actual  length  for  frequencies  less 
than  or  equal  to  the  folding  frequency.  For  data  with  large  holes,  the  effec- 
tive length  gets  shorter  as  the  frequency  increases. 

The  data  spacing  factor  can  be  viewed  as  the  result  of  a  discrete 
integration  over  the  area  where,  but  not  including,  t.  =  t.  to  the  boundary 
t.  =  t.  +  r  .  See  Figure  Bl  for  the  nomenclature.  Some  area  outside  the 
squares  wholly  in  the  integration  area  will  be  included  with  this  approach. 
This  makes  the  integration  area  larger  than  it  actually  is,  leading  to  a  dis- 
torted value  for  the  data  spacing  factor.  Hence  the  following  adjustment  is 
made  in  the  equation  (5). 

F(t. ,t  .)  =   5.8  .    -  A.  .  [9] 

The  A. .  represents  the  area  to  be  subtracted  from  the  rectangular  local  area. 
The  following  diagnostic  criteria  are  needed: 

P  :     If  t..  ^  ^w-i  ^  ^  '  ^^^-^   Pp  i^  inside  the  domain  of   integration: 
otherwise  P  is  outside. 

P  •     If  t.,  .  -  t..  .  +  T  ,  then  P  is  inside  the  domain  of  integration; 
otherwise  both  P  and  P  are  outside. 

P  :      If  t  . .  ,  ^  t . ,  .  +  r  ,  then  P   is  inside  the   domain  of   integration; 

otherwise  both  P„  and  P,  are  outside 
2      4 

P  :      If  t  .    -   ^-1        "•"  ^  '  then  P   is  inside  the  domain   of   integration; 

otherwise  P, ,  P„ ,  P^ ,  and  P^  are  outside. 
12   3       4 

These  criteria  define  the  area  cut  off  the  local  rectangle  by  the  diagonal 
line  t  =  t .  +  T  .  A  point  on  the  boundary  is  defined  to  be  inside  the 
domain.   For  convenience  the  following  notation  is  used: 

d  .     length  cut  off  the  rectangle  in  the  i-direction  along  the  side  meas- 
ured from  point  P 

d       lengt'.  cut  off  in  the  j-direction,  measured  from  point  P  ,  etc. 
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4-1 


(ti_iHi)/2  - 


ti^ 


lVW/2-.(P3> 


4+1 


\ 


T 1 ^-r 


tfi 


Vi  'Vi'V  ^i      ^VVi^ 


V^i^. 


t. 


i+1 


Figure  Bla.  Nomenclature  for  a  rectangular  local  area  associated  with  a 
covariance  ordinate. 
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Figure  Bib.   Example  of  partial  rectangular  local  areas  at  domain  boundaries. 

Source:   Marquardt  and  Acuff  (1985,  Figure  la  and  lb). 

When  t.-t.  >  r  ,  the  areas  of  these  partial  rectangles  does  not   need  to 
be  computeii  since  ?heir  covariance  ordinate  location  (t  ,t  )  lies  outside  the 


56 

domain  of  integration,  where  the  lag  window  D(r)  is  zero.   Where   the  covari- 
ance  lies  inside  the  domain  of  integration,  the  following  cases  could  apply. 


Case  1:   Pp  ■'^^  inside 


A.  .  =  0 


Case  2:   P   is  outside.   P   P   and  P  are  inside.  Then 

d-.  =  t.,  ,-t.,.      >0. 
2j    J(+)     i(-)    m 

A.  .  =  dJ./2. 
Case  3:   P  and  P  are  outside.   P  and  P  are  inside.   Then 


'2j  =  ^j(.)  -  tj(-,  >  0 


^21  =  ^j(.)  -  ^i(-)    -  ^m>«- 

\j  =  ^^i  ^'^2i)^2j/2- 

Case  4:   P  and  P  are  outside.   P  and  P  are  inside.   Then 

d-.  =t.,  ,  -t.,  ,  >0 
21    i(+)    i(-) 


d-.  =  t.,  ,-t.,,-r  >0 
2  J    j(  +  )    i(-)    m 


d.  .  =  t.,  ,-t.,  .-T  >0 
4j    J(+)     i(+)    m 


^j  =  ^'23    ^   '^4j"*2i/2. 


Case  5:   P   P   and  P  are  outside.   P   is  inside.   Then 

d„.  =t.,,-t.,  ,>0 
21     i(+)     i(-) 


^2j  =  ^j(.)  -  tj(-)  >0 
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11 


=  ^j(-l 


^i(- 


>  0 


4j 


^J(^) 


-  t 


i(  +  ) 


>  0. 


A.  . 
IJ 


=  ^2i'^4j  ^  ^^li^^2i» 


2j 


-  d^  .)/2. 


Case  6:  If  the  period  of  observation  begins  before  the  first  time  point  (t  > 
0) ,  or  if  t  <  T,  the  local  areas  at  the  ends  of  the  domain  of 
integration  are  given  by  the  definitions  for  the  endpoints  and  bounded 
by  equation  (6) . 

Case  7:  The  covariance  ordinates  on  the  main  diagonal  where  i  =  j  are  not 
used,  although  some  of  the  area  is  within  the  integration  area.  The 
triangular  region  above  the  main  diagonal  can  be  associated  with  the 
areas  of  the  "above"  and  "right"  rectangles.  The  area  is  split 
between  these  two  regions  by  perpendiculars  from  the  main  diagonal 
connecting  to  the  point  P  of  the  rectangle  above  in  the  integration 
area  and  from  the  main  diagonal  to  the  point  P  of  the  triangle  to  the 
right.  When  the  area  is  unaffected  by  the  restraints  of  the  minimum 
time  increment  for  the  Nyquist  folding  frequency,  the  points  P  and  P 
coincide.  Case  7  applies  only  to  the  points  on  the  main  diagonal  with 
the  associated  rectangles  for  which  j  =  i  +  1.  The  data  spacing  fac- 
tor is  then  modified  by  adding  the  area  of  the  triangles. 


F(t^.t.) 


S.Sj    -  A.  .  + 
1       IJ 


B. 


ij 


when  j-i+1 


[lOj 


There  are  five  possible  subcases. 

7.1 

t.    +  t.  t  •  +  t . 

''   h(.)    <  "2   '   -^"^/^^   S(-)  '     '      2'" 

then  the  point  P  for  the  rectangle  is  constrained  away  from  the 
main  diagonal  and  B.  .  =  0. 


7.2 


If  t.     +  (t  .,  .  -  t .,  .  )/2  <  t.  +  7- 
i(  +  )     j(-t-)     J(-)       1   4<o 


then  the  entire  triangular  area  below  the  rectangle  is  appended  to 

the  rectangle ,  and  B..  =  (t.,  ,  ~  t.,  ,)  /4. 

ij     J(+-)    J(-) 


7.3 


If  t  .,  .  -  (t. ,  ,  -  t. ,  ,)/2  >  t .  -  7- 
J(-)      l(  +  )     i(-)        J    4(0 
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then  the  entire  triangular  area  to  the  left  of  the  rectangle   is 
apporiued  to  the  rectangle,  and  B.^  -  (t     ■  "" ;  /  _  i  ^  /'*- 

7.4 

then  a  trapezoidal  area  below  the  rectangle  is  appended  and 

7.5 

then  a  trapezoidal  area  to  the  left  of  the  rectangle   is  appended 
and 

B.  .  =  (t.,  ,  -  t.  +■  7")  (t.  -  -f-  -  t.  .  ,  ) 
ij     J(-)     J   4(0    J    4(0    i(-) 

To  follow  the  recommendation  of  Marquardt  and  Acuff,  the  spectral  window 
is  evaluated  in  parallel  with  the  estimation  of  the  spectrum.  The  window  has 
the  shape 

n  n         n  n 
Z*((o)  -;p5:EW.  .  =p-ZZ  D(r)  F(t..t  )  cos  27:<o(r).  [12] 

i<j   ^■^      i<j 

When  working  with  evenly  spaced  data.  Z*((o)  is  a  periodic  function  of  w,  the 
sharp  center  lobe  being  repeated  at  all  even  multiples  of  the  Nyquist  fre- 
quency. This  property  leads  to  100  percent  aliasing  of  all  frequencies  above 
the  folding  frequency.  When  dealing  with  irregularly  spaced  data,  Z*(m)  is  an 
almost  periodic  function,  with  the  possibility  of  having  important  side  lobes 
far  from  the  center  lobe  and  within  the  principal  frequency.  Due  to  the 
nature  of  irregularly  spaced  data,  the  spectral  window  has  two  moderate  defi- 
ciencies when  dealing  with  irregular  instead  of  even  data: 

1.  The  center-lobe  is  less  sharp,  which  leads  to  somewhat  higher  variances 
of  the  spectral  estimates. 

2.  Care  must  be  taken  in  interpreting  the  estimated  spectra  due  to  the  pos- 
sibility of  partial  aliasing  if  important  side  lobes  are  present  in  the 
actual  window. 

The  weights  in  the  DQSE  quadratic  form  will  change  with  each  set  of  data  due 
to  the  change  in  spacing.  They  can  be  adjusted  to  sharpen  the  shape  to  yield 
a  good  window.  Characteristics  of  a  good  window  include  a  sharp  center  lobe, 
minimum  area  outside  the  center  lobe,  being  positive  almost  everywhere,  by 
constraint  if  necessary,  and  sometimes  smoothness.  Marquardt  and  Acuff  (1982) 
show  examples  of  the  DQSE  method  without  the  B .  ^,  factor  and  discuss  both  the 
soectrum  and  the  window. 
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APPENDIX  C 
ACTUAL  AND  REGENERATED  DATA 
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APPENDIX  D 
FORTRAN  PROGRAM  TO  REGENERATE  TIME  SERIES 

cccccccccccccccccccccccccccc 
c 

C  THIS  PROGRAM,  DARMA.F,  GENERATES  AN  ARMA(6,5)  PROCESS  USING 

C  PARAMETERS  COMPUTED  BY  THE  MINIMIZATION  OF  THE  -2  LN( LIKELIHOOD) 

C  FUNCTION  OR  ANY  OTHER  TEST  NUMBERS 

C 

C  AUTHOR  -  NAOMI  REGIER 

C  AGRICULTURAL  ENGINEERING 

C  SEATON  HALL 

C  KANSAS  STATE  UNIVERSITY 

C  MS  CLASS  OF  1986 

C 

C  IMPORTANT  NOTE:   THE  ROUTINES,  "RANNP"  AND  "RANN"  ARE  FROM 

C  HARRIS  FORTRAN  AND  MAY  NOT  BE  AVAILABLE  ON  ALL  SYSTEMS. 

C 

C  THE  SUBROUTINE  RANNP  INITIALIZES  THE  ARRAY  FOR  RANN. 

C  ITS  ARGUMENT  IS  THE  SEED,  BY  DEFAULT  =  1. 

C 

C  THE  FUNCTION  RANN  IS  A  REAL  RANDOM  NUMBER  DISTRIBUTED  ABOUT 

C  A  NORMAL  DISTRIBUTION  CURVE.   ITS  FIRST  ARGUMENT  IS  THE  MIDPOINT 

C  OF  THAT  CURVE  AND  THE  SECOND  ARGUMENT  IS  THE  STANDARD  DEVIATION. 

C 

C  Input  is  UNIT  7 

C  Output  is  UNIT  8 

C 

C  VARIABLE  LIST 

C        INPUT  VARIABLES: 

C  RMEAN  -  Desired  mean  of  the  regenerated  numbers 

C  DEV  -  Standard  deviation  of  the  random  error  input 

C  NP  -  Number  of  autoregressive  parameters 

C  NQ  -  Number  of  moving  average  parameters 

C  A  -  Coefficients  for  autoregressive  numbers 

C  B  -  Coefficients  for  moving  average  numbers 

C        OUTPUT  VARIABLES: 

C  A,  B  -  Provide  label  of  coefficients  at  top  of  file 

C  XM  -  Predicted  number  with  a  nonzero  mean 

C        INTERNAL  VARIABLES: 

C  COEF  -  Array  of  random  values,  normally  distributed, 

C  at  time  t-1  to  t-q 

C  X  -  Array  of  past  predicted  values  back  to  time  p 

C  RNEW  -  Random  value  at  time  t 

C  XSUM  -  Amount  due  to  autoregressive  contribution 

C  CSUM  -  Amount  due  to  moving  average  contribution 

C  XNEW  -  Curent  predicted  value 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

REAL  X(6),A(6),B(5),COEF(5) 
C  Initialize  gaussian  random  number  generator 
CALL  RANNP(3.) 
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C  Read  input 

READ(7,*)  RMEAN 
READ(7,*)  DEV 
READ(7,*)  NP,  NQ 
C  Read  in  autoregressive  coefficients 
DO  5  1=1, NP 
5  READ(7,*)  A(I) 
C  Read  in  moving  average  coefficients 
DO  10  J  =  l.NQ 
10  READ(7,*)  B(J) 
C  Write  out  coefficients 
WRITE( 8,2000)  A,B 
2000FORMAT('  ','A:  ',6(F6.3 ,1X) , '  B:  '.5(F6.3 ,1X)) 
C  Get  random  numbers  for  moving  average  noises 
DO  15  I  =  l.NQ 

COEF(I)  =  RANN(0.,DEV) 
15  CONTINUE 
C  Initialize  X  values  with  zeroes 
DO  20  J  =  l.NP 
X(J)  =  0.0 
20  CONTINUE 
C  Recursively  generate  series  of  numbers  from  ARMA(p,q)  model 
DO  100  I  =  NP+1,  600 
RNEW  =  RANN(0.,DEV) 
C  Add  up  autoregressive  contribution  to  present  value 
XSUM  =  0. 
DO  25  K  =  l.NP 
25    XSUM  =  XSUM  +  A(K)  *  X(K) 
C  Add  up  moving  average  contribution  to  present  value 
CSUM  =  0.0 
DO  30  KK  =  1,NQ 
30    CSUM  =  CSUM  +  B(KK)  *  COEF(KK) 
C  Present  value  with  zero  mean 

XNEW  =  XSUM  +  CSUM  +  RNEW 
C  Adjust  mean  to  actual 

XM  =  XNEW  +  RMEAN 
C  Reset  calculated  values  for  new  calculation 
DO  35  J  =  0,  NP-2 
35    X(NP-J)  =  X(NP-J-l) 
X(l)  =  XNEW 
C  Reset  random  numbers  for  next  calculation  of  X 
DO  40  J J  =  0,NQ-2 
40     COEF(NQ-JJ)  =  COEF(NQ-JJ-l) 
COEF(l)  =  RNEW 
C  Discard  first  300  values  to  eliminate  effects  of  initial  zeroes 
IF(I  .GT.  300)  THEN 
WRITE( 8,1000)  XM 
1000       FORMATCIX,  F10.2) 
END  IF 
100  CONTINUE 
STOP 
END 
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ABSTRACT 

Evenly-spaced  time  series  with  missing  values  of  tractor  power  loading 
variations  were  analyzed  by  spectral  density  and  autoregressive-moving  average 
(ARMA)  methods.  The  spectral  density  method  was  unsatisfactory  for  easily 
reproducing  a  series  on  the  dynamometer  and  so  the  coefficients  in  the  ARMA 
model  were  found  by  optimizing  the  value  of  -2  ln( likelihood)  and  selecting 
the  order  from  among  the  top  five  models. 

The  series  for  the  different  tillage  tools  show  some  difference  between 
the  implements.  The  drill  and  drag/harrow  have  fairly  small  variations 
between  one  point  in  time  and  the  next.  The  disk  samples  generally  show  a 
moderately  rougher  pattern,  with  variation  between  the  samples.  The  addition 
of  a  springtooth  behind  the  disk  did  not  produce  any  distinguishable  differ- 
ence. The  chisel  shows  the  roughest  series  while  the  field  cultivator  is  only 
slightly  smoother.  When  NHS  is  applied,  the  power  increments  from  one  value 
to  the  next  are  small  although  the  general  trend  may  vary. 

Appendices  include  an  outline  of  the  basic  theory  for  continuous  and 
discrete  evenly,  and  irregularly-spaced  spectral  density  analysis,  as  well  as 
graphs  of  actual  data  and  the  series  regenerated  from  the  ARMA  models.  The 
algorithm  for  the  ARMA  model  with  missing  values  is  given  in  the  main  part  of 
the  text . 


